Venkat Subramanian

386 Reputation

13 Badges

15 years, 247 days

MaplePrimes Activity


These are replies submitted by

@ChrisRackauckas 

 

Until problems have been benchmarked, the table is not justified. For example, for every row you need to give an example and justify. The previous version of the table and the ranking on Maple seems to be made in a rush doubting the validity of other entries as well. I am not well versed with the other platforms to accept or deny your conclusions. I think, the best improvement I see is the inclusion of updated linear algebra solvers, but it remains to be seen how well it works for largescale problems, stringent tolerances, etc. 

Your table ranks sensitivity analysis, DASPK should be compared and my guess is that it should be excellent. Same with JACOBIAN from Paul Barton.

Regarding your comment on 10,000 or more ODEs, are you talking about well defined system with well defined pattern for Jacobian or matrix or arbitrarily working with a random set of matrix? I agree that Maple is weak for this, but MATLAB switches to sparse solvers for large systems.

I have had bad experience with shooting methods (in particular single shooting) for BVPs. It will be trivial to break any such code for DAE BVPs.

For BVPS, it may be good to include Cash's solvers and Asher's COLSYS, etc. I am not sure on what basis you have given Fair for Maple for BVPs and Good for MATLAB and Mathematica? If it is trivial for you to add a wrapper, add MEBDF family of sovlers from Cash (it can handle index 2 DAEs very well)

You have given Excellent for SUNDIALS for implicitly defined DAE solvers. What is the index here? RADAU handles DAEs better than anything from SUNDIALS (missing a trivial transformation is not a weakness).

 

 

It was a good discussion. 

As you might know IDA came from DASSL/DASKR. Among all the BDF type solvers for index 1 DAEs, DASPK and DASPK3.0 may be the best. Check, http://www.cs.ucsb.edu/~cse/software.html. It has more options for initialization of DAEs.

Also, for sensitivity analysis, DASPK 3.1 was probably the first solver (BDF type). If not the first, it should be pretty damn good.

If stochastic ODEs are to be solved, I would look at StochSS of Petzold before concluding that other software may be the best.

The first successful implementation of sparse solver with BDF type solver was done by Paul Barton at MIT. Check, DAEPACK, and JACOBIAN at https://yoric.mit.edu/software.

In your blog, I have seen claims of SDIRK type methods doing better than other algorithms, I would like to see some examples if possible. 

Thanks

 

is not properly implemented (for nonlinear ODEs). The documentation is wrong. 

Search for LSODE venkat subramanian mapleprimes. The page is taking too long to load. Help files need a lot of work in Maple.

@ as of now, no. but simple strategies can help convert multipoint BVPs to a single domain BVP. in some, cases, it is not possible.

On an another note, what about DAE BVP? As of now, only COLDAE can handle it.

You may be underestimating the linear solver capability of Maple. I used its sparse solver and we have submitted a paper based on higher derivative methods (only for BVPs). Check https://www.mapleprimes.com/posts/208499-A-New-And-Efficient-Code-In-Maple-For.

Let me know if current implementations of BVPs (not DAEs) in Julia or anything else beats its speed (tested with Maple 2017 which will use all the cores if needed for largescale problems or with more nodes).

@ChrisRackauckas 

 

Note that I commend your progress in Julia, and I will update its failure on some of our models (when matrix size becomes too large) soon.

Regarding the table, I think it is an unfair table without doing a fair comparison. You need to give a model or example to see how different platforms perform.

You give MATLAB excellent for complex numbers, I don't think MATLAB can handle nonlinear complex ODEs using BDF type formula. I am not sure why you think Maple won't solve any complex ODE that other platforms can solve. The comparison should be made based on the ability to solve a given model, not based on the implementation of a particular algorithm. Note that dsolve/numeric if properly used can be tricked to solve any RK method (explicit) by giving the Runge kutta coefficient in proper form as long as you give two formula for prediction (to calculate the error)

The main complaint I have with Maple - not integrating sparse solvers into its ODE solvers and not using Newton Raphson in any of its dsolve/numeric solvers.

Regarding events, I think I didn't make it clear. The current documented events page in Maple handles much more than what MATLAB can do. There are much more options builtin, but not documented yet.  Check https://www.mapleprimes.com/questions/125273-Dsolve-Events-How-To-Control-For-A-Sign-Change

For state-dependent delay, I believe Maple is a leader better than current codes of Hairer at least for one difficult problem,

See, https://www.mapleprimes.com/posts/201287-Delay-Differential-Equations--Maple-Is-A-Leader

Please run this example in Julia and your approaches and get back to me with details on cpu time, etc.

I have run IDA, RADAU, etc. For a very small number of ODEs (less than 100 or so), Maple's dsolve numeric will be as fast as anything else for non-stiff and moderately stiff systems. 

 

Based on my experience,

Maple can solve state dependent delay differential equations, Mathematica can't

The table makes me conclude that the author has no clue.

 

BVP solvers are available in Maple as well and marked as unavailable in the table.

Maple's dsolve solver is faster than MATLAB (in compiled form) at least till 200 or more coupled ODEs.

Event options are better in Maple compared to matlab or anything else in that table, but poorly documented.

Maple doesn't have a true DAE solver, but can handle simple index 2 and index 3 systems.

The comment on wrapper for lsode is wrong. Maple has a wrong implementation of LSODE (linearized Newton Raphson).

Maple can be used to solve ODEs with arbitrary precision, also ODEs in the complex domain. It looks the blogger hasn't used a recent version of Maple.

 

@rlopez 

I have not been to install and use classic worksheet in windows10 computers and no one at Maplesoft has helped me in fixing the same. The software gets installed, but when I open it, I can't save any file and the window disappears. 

This is for various Maple versions, from Maple 14 to Maple 2017 running in all possible compatible modes.  This means that Maplesoft probably didn't test and validate the robustness of classic worksheets in windows 10 yet.

 

My group is still using windows 7 for just one reason. The need to use classic worksheet Maple.

 

@Preben Alsholm 

Try epsilon = 10, in the example I gave. HDM works, but maxmesh = 512 won't work. Newton iteration is failing. So, continuation or use of an approximation solution or higher initmesh/maxmesh will be needed.

Regarding a problem with multiple solutions, HDM can handle it by calling HDM (instead of HDMAdapt) with different initial guesses. Please check the uploaded pdf listed for this.

PS, it is not my intention to make this thread a discussion on when dsolve/numeric fails for BVPs and how to fix it. (I do believe that this should be done as a separate thread with detailed examples as current help files do not justify the capability of dsolve/numeric). dsolve/numeric has all the bells and whistles, but HDM, as written, is just a functional BVP solver comparable to Fortran, MATLAB or C routines available for BVPs in terms of ease of use. That said, because of the symbolic capability of Maple, higher derivatives are easily found, and combined with the Sparse solver, HDM has turned out to be competitive. For large scale problems, Maple's sparse solver in recent versions of Maple even uses all the processors in parallel to gain additional speed.

 

@Preben Alsholm 

try epsilon = 8, at default settings, dsolve/numeric fails. Of course, it is possible to increase maxmesh, introduce continuation, etc. But at the default setting, HDM is better as it has 10th order accuracy. Also, see 35 examples of Cash posted in my website, dsolve/numeric failed for many of the cases without continuation or changing the settings. Some problems couldn't be solved even after making the changes. (Try very low values of the perturbation parameters for the problems given by  Cash).


 

Reset the program to clear the memory from previous execution command.

restart:

 

Read the txt file which contains the HDM solver for BVPs.

read("HDM.txt");

 

Declare the precision for the entire Maple® sheet.

Digits:=15;

Digits := 15

(1)

 

Enter the first-order ODEs into EqODEs list.

EqODEs:=[diff(y1(x),x)=y2(x),diff(y2(x),x)=epsilon*sinh(epsilon*y1(x))];

EqODEs := [diff(y1(x), x) = y2(x), diff(y2(x), x) = epsilon*sinh(epsilon*y1(x))]

(2)

 

Define the left boundary condition (bc1), and the right boundary condition (bc2). One should collect all the terms in one side.

bc1:=evalf([y1(x)]);

bc1 := [y1(x)]

(3)

bc2:=evalf([y1(x)-1]);

bc2 := [y1(x)-1.]

(4)

 

Define the range (bc1 to bc2) of this BVP.

Range:=[0.,1.];

Range := [0., 1.]

(5)

 

List any known parameters in the list.

pars:=[epsilon=8];

pars := [epsilon = 8]

(6)

 

List any unknown parameters in the list. When there is no unknown parameter, use [ ].

unknownpars:=[];

unknownpars := []

(7)

 

Define the initial derivative in nder (default is 5 for 10th order) and the number of the nodes in nele (default is 10 and distributed evenly across the range provided by the user). The code adapts to increase the order. For many problems, 10th order method with 10 elements are sufficient.

nder:=5;nele:=10;

nder := 5

nele := 10

(8)

 

Define the absolute and relative tolerance for the local error. The error calculation is done based on the norm of both the 9th and 10th order simulation results.

atol:=1e-6;rtol:=atol/100;

atol := 0.1e-5

rtol := 0.100000000000000e-7

(9)

 

Call HDMadapt procedure, input all the information entered above and save the solution in sol. HDMadapt procedure does not need the initial guess for the mesh.

sol:= HDMadapt(EqODEs,bc1,bc2,pars,unknownpars,nder,nele,Range,atol,rtol):

"HDM fail, increase points to...", 20

"HDM fail, increase points to...", 40

"HDM fail, increase points to...", 80

"case2 - insert/remove points based on 2n-1 error to...", 48

"case1 - double the nodes to...", 96

(10)

 

Present some details of the solution.

sol[4]; # final derivative

5

(11)

sol[5]; # Maximum local RMSE

0.430441018168833e-7

(12)

 

Store the dimension of the solution (after adjusting the mesh) to NN.

NN:=nops(sol[3])+1;

NN := 97

(13)

 

Plot the interested variable (the ath ODE variable will be sol[1][i+NN*(a-1)] )

node:=nops(EqODEs);
odevars:=select(type,map(op,map(lhs,EqODEs)),'function');

node := 2

odevars := [y1(x), y2(x)]

(14)

xx:=Vector(NN):

xx[1]:=Range[1]:

for i from 1 to nops(sol[3]) do xx[i+1]:=xx[i]+sol[3][i]: od:

for j from 1 to node do
  plot([seq([xx[i],rhs(sol[1][i+NN*(j-1)])],i=1..NN)],axes=boxed,labels=[x,odevars[j]],style=point);
end do;

 

sol1:=dsolve({op(subs(pars,EqODEs)),y1(0)=0,y1(1)=1},type=numeric):

Error, (in dsolve/numeric/bvp) uanble to achieve continuous solution with requested accuracy of .1e-5 with maximum 128 point mesh (was able to get .21e-3), consider increasing `maxmesh` or using larger `abserr`

 

sol1(1);

sol1(1)

(15)

 

 

 


 

Download Example_3_Troesch.mws

@JohnS 

Thanks for the reply and support. Please note that this code fails less frequently compared to dsolve/numeric(based on our testing). In addition, mixed boundary conditions can be easily addressed. There are many ways. One way is to add a dummy variable for every mixed boundary condition as done in Maple's dsolve/numeric.

For example if you have

dy1/dx = -y1^2

dy2/dx = -0.1y2^2

with bcs y1(0)-1=0

y2(0)-y1(1)=0

then you add dy3/dx =0 and rewrite as

dy1/dx = -y1^2

dy2/dx =-0.1y2^2

dy3/dx = 0

y1(0)-1=0

y2(0)-y3(0)=0

y1(1)-y3(1)=0

 

But our tests indicate that this may not be the best way in particular for a problem with a large number of mixed bcs as y3 is solved at every node. So, once the code is optimized it will be released. Right now you can use the approach mentioned in this thread.

CAM-D-17-01905.pdf

HDM.txt (updated 08/26/17, a small typo was causing issues for some problems).

 

(Ideally, these files should have been included with the main post, but I couldn't do it).

if you like the classic worksheet, say good bye to it in windows 10. I have not been able to create, save or do anything in the classic worksheet (from maple 14 to 2017) in windows 10.

The weird thing is if I upgrade windows 7 laptops to 10, some versions of classic worksheet work. For new laptops that come pre-installed with windows 10, I have had no luck with having classic worksheet stay open for even 1 minute.

 

@ code as given gives meaningful results till time = 0.1 with the method chosen. Default method fails even till t= 0.1.

@Markiyan Hirnyk 

see method option for pdsolve/numeric, and one of that might (should) work for the separated PDE. Some attempts are given below.
 

``

restart;

sys:={diff(u(x, t), t)-diff(v(x, t), x)+u(x, t)+v(x, t) = (1+t)*x+(x-1)*t^2, diff(v(x, t), t)-diff(u(x, t), x)+u(x, t)+v(x, t) = (1+t)*x*t+(2*x-1)*t};
##
icbs:= {u(0, t) = 0, u(x, 0) = 0, v(0, t) = 0, v(x, 0) = 0};
##
pde1:=eval(sys[1]-sys[2],u(x,t)=v(x,t)+w(x,t));
solw:=pdsolve(pde1,{w(x,0)=0,w(0,t)=0},numeric,time=t,timestep=0.01,spacestep=0.01,range=0..1);
solw:-plot3d(x = 0 .. 1, t = 0 .. 1); #No problem
pde2:=eval(sys[1]+sys[2],u(x,t)=-v(x,t)+z(x,t));
solz:=pdsolve(pde2,{z(x,0)=0,z(0,t)=0},numeric,time=t,range=0..1,timestep=0.01,spacestep=0.01);
solz:-plot3d( x = 0 ..1, t = 0 .. 1); #Problem

sys := {diff(u(x, t), t)-(diff(v(x, t), x))+u(x, t)+v(x, t) = (1+t)*x+(x-1)*t^2, diff(v(x, t), t)-(diff(u(x, t), x))+u(x, t)+v(x, t) = (1+t)*x*t+(2*x-1)*t}

icbs := {u(0, t) = 0, u(x, 0) = 0, v(0, t) = 0, v(x, 0) = 0}
pde1 := diff(w(x, t), t)+diff(w(x, t), x) = (1+t)*x+(x-1)*t^2-(1+t)*x*t-(2*x-1)*t

solw := _m220916848

 

 

pde2 := diff(z(x, t), t)+2*z(x, t)-(diff(z(x, t), x)) = (1+t)*x+(x-1)*t^2+(1+t)*x*t+(2*x-1)*t

solz := _m219596616

Error, (in pdsolve/numeric/plot3d) unable to compute solution for t>HFloat(0.26000000000000006):
solution becomes undefined, problem may be ill posed or method may be ill suited to solution

 

 

solz:=pdsolve(pde2,{z(x,0)=0,z(0,t)=0},numeric,time=t,range=0..1,method = ForwardTime1Space[backward], timestep = 0.00001,compile=true);
solz:-plot3d( x = 0 ..1, t = 0 .. 0.1,axes=boxed);

solz := _m195442264

 

 

 

NULL


 

Download pdehyperbolic.mw

 

5 6 7 8 9 10 11 Last Page 7 of 17