Advanced Search

Venkat Subramanian

431 Reputation

13 Badges

15 years, 265 days
Contact 119368

MaplePrimes Activity


These are replies submitted by

Definitely wrong...

Dr. Venkat Subramanian 431
September 26 2017

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.

 

To me it appears that way...

Dr. Venkat Subramanian 431
August 28 2017
0 0

@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.

 

Problem dependent...

Dr. Venkat Subramanian 431
August 26 2017

@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...

Dr. Venkat Subramanian 431
August 25 2017

@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

It is straightforward...

Dr. Venkat Subramanian 431
August 24 2017

@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.

Not able to upload the pdf and text file...

Dr. Venkat Subramanian 431
August 23 2017

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).

classic won't work...

Dr. Venkat Subramanian 431
July 21 2017
0 2

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.

 

tweak node spacing, time step and method...

Dr. Venkat Subramanian 431
April 30 2017

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

perhaps you might want to read What I wr...

Dr. Venkat Subramanian 431
April 30 2017

@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

 

no bug, just not implemented for multipl...

Dr. Venkat Subramanian 431
April 28 2017

@Markiyan Hirnyk 

 

Maple provides different numerical methods for a single PDE (including methods that can handle hyperbolid PDEs in which convection dominates). Your example is a convective system and Maple's current implementation (if I am not wrong) is central difference in x and simple forward scheme in time. This is not stable for convective PDEs and hence you have this error.

That said, I have seen Maple failing even a single diffusion type PDE (linear), I will report that later. In particular, my concern was when it gave wrong numerical results (I can take a failure as in this case, but giving wrong results is dangerous).

Thanks...

Dr. Venkat Subramanian 431
December 08 2016
0 0

@Joe Riel 

Is it possible to delete a file and add a modifed file (with print inserts) to the same library with the same name?

 

can you attach a working code?...

Dr. Venkat Subramanian 431
December 08 2016
0 0

@ 

thanks

 

is this in regular or classic worksheet ...

Dr. Venkat Subramanian 431
December 08 2016
0 0

@Preben Alsholm 

It is not working in classic worksheet. thanks

 

Note that eigenvalues may not be always ...

Dr. Venkat Subramanian 431
November 28 2016
0 0

@Preben Alsholm 

 

I suggest direct dsolve or MatrixExponential(A,t) as the general solution. A can be singular, might have complex eigenvalues, repeated eigenvalues, etc. 

Only when A has distinct eigen values, one can use the lamba,P approach

No,...

Dr. Venkat Subramanian 431
November 26 2016
0 2

@Giulianot 

Not all DAEs can b converted to da/dt = something, db/dt = something. That is the issue.

First 6 7 8 9 10 11 12 Last Page 8 of 18