@vv   AFAIK, numeric dsolve never automatically chooses the method other than distinguishing among IVPs, BVPs, and DAE IVPs. Do you have evidence otherwise?

If it's an IVP that hasn't explicity been declared stiff, and no method is declared, then the method will be rkf45.

@ajfriedlan   If this rkf45 solution runs too slow for you (due to the high accuracy requirement), it may be worthwhile specifying another method, such as dverk78.

@Subhan331 The current example is a nonlinear function. Try changing the line where the function is defined. 

And why would you want to directly specify the number of iterations? I could see wanting to limit the maximum number of iterations to prevent infinite loops, but what's the point of extra iterations after it's converged to your desired error tolerance? 

@Subhan331 Is there something wrong with the code that you posted above? I could suggest some minor improvements, but it is a correct implementation of Newton's method.

@mehdibaghaee The ^+ is a transpose operator which I guess doesn't work in your version of Maple or in 2D Input. You can replace it with ^%T

plot([seq(PPP[[1,k],..]^%T, k= 1..rhs(rtable_dims(PPP)[1]))], color= black);

@mehdibaghaee Redundant or tedious commands are almost never needed in Maple. Your sequence of commands above can be reduced to the single line

plot([seq(PPP[[1,k],..]^+, k= 2..rhs(rtable_dims(PPP)[1]))], color= black);

Note that there's no need to specify with(plots), the number of plots (24 in this case), or the number of points per plot (g in this case).

@phbmc What range of values of b are you interested in? For b < 10^5 (approximately), the "brute force" method works well, and I'll post some code soon. That just means trying all possible a1 and checking if the corresponding a2 is an integer in the correct range.

@Subhan331 The code says that the loop ends when the absolute error is less than tol. It just so happens in this case that that happens on the sixth iteration. If you want more iterations, use a smaller value of tol.

As -log[10](tol) gets close to Digits, you will also need to increase Digits to continue to have meaningful iterations. Newton's method tends to double the number of correct Digits on each iteration once it's zeroed in on a target.

What do you want to achieve by placing an index, such as [maximum], on a non-atomic algebraic expression, such as (-k*a)? Maple does allow this, but I find it unwieldy, unless it's meant strictly for display purposes.

@Kitonum Note that the boundary of the region is in general the complicated and disconnected set {(x,y) | f(x,y) = c }. A general treatment of cases of polynomial f might be possible with solve(..., parametric). I don't know how far beyond polynomial the capabilities of that command has been extended over the past few years.

@Preben Alsholm Thank you for your comments and your totally nontrivial digging through the dsolve code.

You wrote:

• Maybe I'm misunderstanding your intention, but the localized version still knows about the parameter 'a' and that can still be set:

What I mean by A is an instantiation of X is that X is some structure with a free variable x (also known in mathematical language as a "parameter" (but "parameter" means some else in computer science)), and A is essentially a copy of X except that x has been given a specific value, a. The following two further things are irrelevant to my notion of an instantiation, but I think that you are trying to ascribe them to it:

1. whether or not A can be further manipulated to change the parameter again;
2. whether or not A retains enough of the original information of X that X could in theory be reconstructed (aka "reverse engineered") from A.

By the instantiation being localized, I mean that A has been saved or stored in such a way that if its parent X is instantiated again with a different value of x, this will not affect A.

These terms may have other definitions, and I may be using them incorrectly due to my having come up with the concepts independently (and, if so, I'd appreciate a correction).

Acer: Did you understand my previous usages to have the defintions above?

• I thought that you wanted a procedure like the one you get when using dsolve/numeric for a concrete value of 'a'.

No. If that were feasible, then indeed it would serve my purpose, but I doubt that it's feasible, and I never considered it.

• In fact Sol and res above are identical.

Yes, and that seems irrelevant to me. What's relevant is whether you can use the parent, Sol, to create res1 and res2 such that res1 and res2 exist in the same process's RAM at the same time and such that res1 <> res2 and furthermore they are not just superfically different (like Matrices with the same entries), but essentially different in that they act with different parameter values. You discovered a way to do that by using .m files. My technique does the same thing with %m formats, entirely in RAM. Yes, those procedures are quite long---this method does incur some significant overhead. Nonetheless, there is a substantial time savings.

• Furthermore, if you use the "naive" approach, but add option remember it seems that you get good results for the plot3d problem.

It's only happenstance that your code works (produces the correct plot) because it's a sheer coincidence that plot3d happens to access the GRID coordinates in the order that minimizes the number the changes of the parameters. If you change the plot to

plot3d(Y(A,X), X= -2..2, A= -2..2, ...),

you'll see that your code gets it wrong while mine doesn't. In particular, yours acts as if the parameter remains unchanged after the final A value is first used.

Yes, instantiations can be localized, completely in RAM, using only the half-line of code given in the title and a remember table! This gives a brief solution to the timing anomaly that is the subject of this Post. Specifically, we can guarantee that dsolve's parameters will be set exactly once for each numeric valuation of the parameters, regardless of the order that plot3d (or any other code) chooses to evaluate the (independent variable, parameter) pairs.

restart: 
Sol:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):

#Naive method:
Y:= (A,X)-> 
   if [A,X]::[numeric$2] then Sol('parameters'= [A]); eval(:-y(:-x), Sol(X))
   else 'procname'(args)
   fi
:
gc():
P1:= CodeTools:-Usage(plot3d(Y(A,X), A= -2..2, X= -2..2, grid= [99$2])):

#Save data array to verify that the other method gives the same plot.
M1:= op([1,3], P1):
save M1, "/Users/carlj/desktop/M1.m":

memory used=325.30MiB, alloc change=4.00MiB, cpu time=2.84s, real time=2.84s, gc time=140.62ms


restart:
Sol:= dsolve({diff(y(x),x)=a*y(x), y(0)=1}, numeric, parameters= [a]):

#"Carl's Magic .m" method:
Localize:= proc(Sol, params::list(complexcons))
option remember;
   Sol('parameters'= params);
   sscanf(sprintf("%m", eval(Sol)), "%m")[]:
end proc:
Y:= (A,X)-> 
   `if`([A,X]::[numeric$2], eval(:-y(:-x), Localize(Sol, [A])(X)), 'procname'(args)):
gc():
P2:= CodeTools:-Usage(plot3d(Y(A,X), A= -2..2, X= -2..2, grid= [99$2])):

#Compare with saved plot data:
M2:= op([1,3], P2):
read "/Users/carlj/desktop/M1.m":

`Plot difference` = LinearAlgebra:-Norm(Matrix(M1-M2));

memory used=90.03MiB, alloc change=32.00MiB, cpu time=688.00ms, real time=699.00ms, gc time=31.25ms
                      Plot difference = 0.

Thank you for this valuable Post. I wonder if there's a way to localize dsolve's returned solution procedure(s) instantiated at particular parameter values. I tried:

Sol:= dsolve(..., numeric, parameters= [...]);
Sol(parameters= [case 1]);
Sol1:= copy(eval(Sol));
Sol(parameters= [case 2]);
Sol2:= copy(eval(Sol));

But, alas, Sol1 now acts as if the case 2 parameters were in effect. I guess this means that the instantiated parameter values are global?

There are two characters appearing in the statement of Problem B1 that have been transcribed incorrectly. The first looks like numeral "2", and I'm sure that it's supposed to be the "is an element of" symbol. The second looks like "y", but I can't figure out what it's supposed to be.

Here is an example of the same thing done as what Mac Dude called a "module factory", i.e., a procedure that returns a module that implements some object (meant in the sense of Object-Oriented Programming or OOP):

CTS:= proc(Init::integer:= 0)
   module()
   local Counter::integer;
   export
      Inc:= proc(inc::integer:= 1) Counter:= Counter+inc end proc,
      Reset:= proc(init::integer:= Init) Counter:= init end proc,
      Query:= ()-> Counter
   ;
      Counter:= Init
   end module
end proc:

Sammy:= CTS(7):
Sammy:-Inc();
Sammy:-Query();
Sammy:-Reset();
                               8
                               8
                               7

Note that the code in significantly shorter, mostly because the module no longer needs to manage the names of the individual counters. The New has been replaced by the outer procedure itself. Note that the module contains one statement.

In terms of its behavior, this is very close to an object module, like I mentioned earlier. However, the syntax of an object module is significantly different.

@Rohith fsolve is a command that returns a purely numeric complex float solution to a system of equations that has the same number of equations as variables. The system of equations (or single equation) may be specified symbolically (as is the case here), or may be specified by a numeric procedure. When I referred to symbolic computation systems, I was referring to the ability to take a symbolic equation and use a one-line template to construct a procedure that numerically solves the equation. That process is symbolic, even though the final product is numeric.