## Numerical solution of challenging ODE...

Hello everybody.

I'm trying to obtain the numerical solution of a differential equation. Unfortunately, this prove to be quite challenging. I was able to obtain a rough solution using mathematica, but nothing more. The function is strictly increasing (for sure).

Any help is really REALLY appreciated, thanks!

 (1)

 (2)

 (3)

## Delay differential equations - Maple is a leader

Maple's dsolve numeric can solve delay ODEs and DAEs as of Maple 18. However, if I am not wrong, it cannot solve delay equations with a time dependent history. In this post I show two examples.

Example 1:

y1(t) and y2(t) with time dependent history. Use of piecewise helps this problem to be solved efficiently. Hopefully Maple will add history soon in its capability.

Example 2:

This is a very a complicated stiff problem from immunology. As of now, I believe only Maple can solve this (other than RADAR5 from Prof. Hairer). Details and plots are posted in the attached code.

Let me know if any one has a delay problem that needs to be solved. I have tested many delay problems in Maple (they work fine). The attached examples required addtional tweaking, hence the post.

I want to take this opportunity to congratulate and thank Maple's dsolve numeric/delay solvers for their fantastic job. Maple is world leader not because of example1, but because of its ability to solve example 2.

 > restart;

This code is written by Dayaram Sonawane and Venkat R. Subramnian, University of Washington. You will need Maple 18 or later for this. For those who are wanting to solve these problems in earlier versions, I can help them by offering a procedure based approach (less efficient).

Example1 The first example solved is a state dependent delay problem (http://www.mathworks.com/help/matlab/math/state-dependent-delay-problem.html).

 > eq1:= diff(y1(t),t)=y2(t);
 (1)
 > eq2:=diff(y2(t),t)=-y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t));
 (2)

Both y1(t) and y2(t) have time dependent history (y1(t)=log(t) and y2(t)=1/t, t<-0.1). If I am not mistaken one cannot solve this directly using Maple's dsolve numeric command. However, a simple trick can be used to redefine the equations for y1(t) and y2(t) as below

 > eq3:=diff(y1(t),t)=piecewise(t<=0.1,1/t,y2(t));
 (3)
 > eq4:=diff(y2(t),t)=piecewise(t<=0.1,-1/t^2,-y2(exp(1-y2(t)))*y2(t)^2*exp(1-y2(t)));
 (4)

The problem is solved from a small number close to t = 0 (1e-4) to make Maple's dsolve numeric remember the history till t = 0.1

 > epsilon:=1e-4;
 (5)
 > sol:=dsolve({eq3,eq4,y1(epsilon)=log(epsilon),y2(epsilon)=1/epsilon},type=numeric,delaymax=5):
 > with(plots):
 > odeplot(sol,[t,y1(t)],0.1..5,thickness=3,axes=boxed);
 > odeplot(sol,[t,y2(t)],0.1..5,thickness=3,axes=boxed);
 > sol(5.0);log(5.0);1/5.0;
 (6)

Tweaking the tolerances and epsilon will get the solution even more closer to the expected answers.

 >
 >

Example 2

The next problem discussed is very stiff, complicated and as of today, according Professor Hairer (one of the world's leading authority in numerical solutions of ODEs, DAEs), cannot be solved by any other code other than his RADAR (5th order implicit Runge Kutta modified for delay equations, Guglielmi N. and Hairer E. (2001) Implementing Radau IIa methods for stiff delay differential equations. Computing 67:1-12). This problem requires very stringent tolerances. For more information read, http://www.scholarpedia.org/article/Stiff_delay_equations. I can safely say that Maple can boast that it can solve this delay differential equation by using a switch function (instead of Heaviside/picecewise function). Code is attached below and results are compared with the output from RADAR code.  Note that dsolve/numeric is probably taking more time steps compared to RADAR, but the fact that Maple's dsolve numeric solved this model (which cannot be solved in Mathematica or MATLAB[needs confirmation for MATLAB]) should make Maple's code writers proud. It is very likely that we will be trying to submit an educational/research article on this topic/example soon to a journal. For some weird reasons, stiff=true gives slightly inaccurate results.

 > restart:
 >
 (7)
 (8)
 > eq[1]:=diff(y[1](t),t)=-r*y[1](t)*y[2](t)-s*y[1](t)*y[4](t);
 (9)
 > eq[2]:=diff(y[2](t),t)=-r*y[1](t)*y[2](t)+alpha*r*y[1](y[5](t))*y[2](y[5](t))*H1;#Heaviside(t-35);
 (10)
 > eq[3]:=diff(y[3](t),t)=r*y[1](t)*y[2](t);
 (11)
 > eq[4]:=diff(y[4](t),t)=-s*y[1](t)*y[4](t)-gamma1*y[4](t)+beta*r*y[1](y[6](t))*y[2](y[6](t))*H2;#Heaviside(t-197);
 (12)
 > eq[5]:=diff(y[5](t),t)=H1*(y[1](t)*y[2](t)+y[3](t))/(y[1](y[5](t))*y[2](y[5](t))+y[3](y[5](t)));#eq[7]:=y[7](t)=HH(t);
 (13)
 > eq[6]:=diff(y[6](t),t)=H2*(10.^(-12)*0+y[2](t)+y[3](t))/(10.^(-12)*0+y[2](y[6](t))+y[3](y[6](t)));
 (14)
 > H1:=1/2+1/2*tanh(100*(t-35));H2:=1/2+1/2*tanh(100*(t-197));
 (15)
 > alpha:=1.8;beta:=20.;gamma1:=0.002;r:=5.*10^4;s:=10.^5;
 (16)
 > seq(eq[i],i=1..6);
 (17)
 > ics:=y[1](0)=5.*10^(-6),y[2](0)=10.^(-15),y[3](0)=0,y[4](0)=0,y[5](0)=1e-40,y[6](0)=1e-20;
 (18)
 > #infolevel[all]:=10;
 > sol:=dsolve({seq(eq[i],i=1..6),ics},type=numeric,delaymax=300,initstep=1e-6,abserr=[1e-21,1e-21,1e-21,1e-21,1e-9,1e-9],[y[1](t),y[2](t),y[3](t),y[4](t),y[5](t),y[6](t)],relerr=1e-9,maxstep=10,optimize=false,compile=true,maxfun=0):

note that compile = true was used for efficiency

 > t11:=time():sol(300);time()-t11;
 (19)
 > with(plots):
 (20)
 (21)

Values at t = 300 match with expected results.

 > p[1]:=odeplot(sol,[t,log(y[1](t))/log(10)],0..300,axes=boxed,thickness=3):
 > display({pr[1],p[1]});
 > p[2]:=odeplot(sol,[t,log(y[2](t))/log(10)],0..300,axes=boxed,thickness=3,numpoints=1000):
 > display({pr[2],p[2]});
 >
 > p[3]:=odeplot(sol,[t,log(y[3](t))/log(10)],0..300,axes=boxed,thickness=3):
 > display({pr[3],p[3]});
 > p[4]:=odeplot(sol,[t,log(y[4](t))/log(10)],197..300,axes=boxed,thickness=3,view=[197..300,-9..-5]):
 > display({pr[4],p[4]});
 > p[5]:=odeplot(sol,[t,y[5](t)],0..300,axes=boxed,thickness=3):
 > display({pr[5],p[5]});
 > p[6]:=odeplot(sol,[t,y[6](t)],0..300,axes=boxed,thickness=3):
 > display({pr[6],p[6]});

## Tricking Maple to solve Differential Algebraic...

Solving DAEs in Maple

As I had mentioned in many posts, Maple cannot solve nonlinear DAEs. As of today (Maple 2015), given a system of index 1 DAE dy/dt = f(y,z); 0 = g(y,z), Maple “extends” g(y,z) to get dz/dt = g1(y,z). So, any given index 1 DAE is converted to a system of ODEs dy/dt = f, dz/dt = g1 with the constraint g = 0, before it solves. This is true for all the solvers in Maple despite the wrong claims in the help files. It is also true for MEBDFI, the FORTRAN implementation of which actually solves index 2 DAEs directly. In addition, the initial condition for the algebraic variable has to be consistent. The problem with using fsolve is that one cannot specify tolerance. Often times one has to solve DAEs at lower tolerances (1e-3 or 1e-4), not 1e-15. In addition, one cannot use compile =true for high efficiency.

The approach in Maple fails for many DAEs of practical interest. In addition, this can give wrong answers. For example, dy/dt = z, y^2+z^2-1=0, with y(0)=0 and z(0)=1 has a meaningful solution only from 0 to Pi/2, Maple’s dsolve/numeric will convert this to dy/dt = z and dz/dt = -y and integrate in time for ever. This is wrong as at t = Pi/2, the system becomes index 2 DAE and there is more than 1 acceptable solution.

We just recently got our paper accepted that helps Maple's dsolve and many ODE solvers in other languages handle DAEs directly. The approach is rather simple, the index 1 DAE is perturbed as dy/dt = f.S ; -mu.diff(g,t) = g. The switch function is a tanh function which is zero till t = tinit (initialization time). Mu is chosen to be a small number. In addition, lhs of the perturbed equation is simplified using approximate initial guesses as Maple cannot handle non-constant Mass matrix. The paper linked below gives below more details.

Next, I discuss different examples (code attached).

Example 1: Simple DAE (one ODE and one AE), the proposed approach helps solve this system efficiently without knowing the exact initial condition.

Hint: The code is given with a semicolon at the end of the most of the statements for educational purposes for new users. Using semicolon after dsolve numeric in classic worksheet crashes the code (Maplesoft folks couldn’t fix this despite my request).

Example 2:

This is a nickel battery model (one ODE and one AE). This fails with many solvers unless exact IC is given. The proposed approach works well. In particular, stiff=true, implicit=true option is found to be very efficient. The code given in example 1 can be used to solve example 2 or any other DAEs by just entering ODEs, AEs, ICs in proper places.

Example 3:

This is a nonlinear implicit ODE (posted in Mapleprimes earlier by joha, (http://www.mapleprimes.com/questions/203096-Solving-Nonlinear-ODE#answer211682 ). This can be converted to index 1 DAE and solved using the proposed approach.

Example 4:

This example was posted earlier by me in Mapleprimes (http://www.mapleprimes.com/posts/149877-ODEs-PDEs-And-Special-Functions) . Don’t try to solve this in Maple using its dsolve numeric solver for N greater than 5 directly. The proposed approach can handle this well. Tuning the perturbation parameters and using compile =true will help in solving this system more efficiently.

Finally example 1 is solved for different perturbation parameters to show how one can arrive at convergence. Perhaps more sophisticated users and Maplesoft folks can enable this as automatically tuned parameters in dsolve/numeric.

Note:

This post should not be viewed as just being critical on Maple. I have wasted a lot of time assuming that Maple does what it claims in its help file. This post is to bring awareness about the difficulty in dealing with DAEs. It is perfectly fine to not have a DAE solver, but when literature or standard methods are cited/claimed, they should be implemented in proper form. I will forever remain a loyal Maple user as it has enabled me to learn different topics efficiently. Note that Maplesim can solve DAEs, but it is a pain to create a Maplesim model/project just for solving a DAE. Also, events is a pain with Maplesim. The reference to Mapleprimes links are missing in the article, but will be added before the final printing stage. The ability of Maple to find analytical Jacobian helps in developing many robust ODE and DAE solvers and I hope to post my own approaches that will solve more complicated systems.

I strongly encourage testing of the proposed approach and implementation for various educational/research purposes by various users. The proposed approach enables solving lithium-ion and other battery/electrochemical storage models accurately in a robust manner. A disclosure has been filed with the University of Washington to apply for a provisional patent for battery models and Battery Management System for transportation, storage and other applications because of the current commercial interest in this topic (for batteries). In particular, use of this single step avoids intialization issues/(no need to initialize separately) for parameter estimation, state estimation or optimal control of battery models.

Appendix A

Maple code for Examples 1-4 from "Extending Explicit and Linealry Implicit ODE Solvers for Index-1 DAEs", M. T. Lawder,

V. Ramadesigan, B. Suthar and V. R. Subramanian, Computers and Chemical Engineering, in press (2015).

Use y1, y2, etc. for all differential variables and z1, z2, etc. for all algebraic variables

Example 1

 > restart;
 > with(plots):

Enter all ODEs in eqode

 > eqode:=[diff(y1(t),t)=-y1(t)^2+z1(t)];
 (1)

Enter all AEs in eqae

 > eqae:=[cos(y1(t))-z1(t)^0.5=0];
 (2)

Enter all initial conditions for differential variables in icodes

 > icodes:=[y1(0)=0.25];
 (3)

Enter all intial conditions for algebraic variables in icaes

 > icaes:=[z1(0)=0.8];
 (4)

Enter parameters for perturbation value (mu), switch function (q and tint), and runtime (tf)

 > pars:=[mu=0.1,q=1000,tint=1,tf=5];
 (5)

Choose solving method (1 for explicit, 2 for implicit)

 > Xexplicit:=2;
 (6)

Standard solver requires IC z(0)=0.938791 or else it will fail

 > solx:=dsolve({eqode[1],eqae[1],icodes[1],icaes[1]},numeric):
 Error, (in dsolve/numeric/DAE/checkconstraints) the initial conditions do not satisfy the algebraic constraints   error = .745e-1, tolerance = .559e-6, constraint = cos(y1(t))-z1(t)^.5000000000000000000000
 > ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));
 (7)
 > NODE:=nops(eqode);NAE:=nops(eqae);
 (8)
 > for XX from 1 to NODE do EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff; end do;
 (9)
 > for XX from 1 to NAE do EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do;
 (10)
 >
 > Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};
 (11)
 > Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};
 (12)
 > icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);
 (13)
 > for j from 1 to NAE do
 > EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);
 > end do:
 > Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};
 (14)
 > if Xexplicit=1 then
 > sol:=dsolve(Sys,numeric):
 > else
 > sol:=dsolve(Sys,numeric,stiff=true,implicit=true): end if:
 >
 > for XX from 1 to NODE do a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red); end do:
 > for XX from NODE+1 to NODE+NAE do a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue); end do:
 > display(seq(a||x,x=1..NODE+NAE),axes=boxed);

End Example 1

 >

Example 2

 > restart;
 > with(plots):
 > eq1:=diff(y1(t),t)=j1*W/F/rho/V;
 (15)
 > eq2:=j1+j2=iapp;
 (16)
 > j1:=io1*(2*(1-y1(t))*exp((0.5*F/R/T)*(z1(t)-phi1))-2*y1(t)*exp((-0.5*F/R/T)*(z1(t)-phi1)));
 (17)
 > j2:=io2*(exp((F/R/T)*(z1(t)-phi2))-exp((-F/R/T)*(z1(t)-phi2)));
 (18)
 > params:={F=96487,R=8.314,T=298.15,phi1=0.420,phi2=0.303,W=92.7,V=1e-5,io1=1e-4,io2=1e-10,iapp=1e-5,rho=3.4};
 (19)
 > eqode:=[subs(params,eq1)];
 (20)
 > eqae:=[subs(params,eq2)];
 (21)
 > icodes:=[y1(0)=0.05];
 (22)
 > icaes:=[z1(0)=0.7];
 (23)
 > solx:=dsolve({eqode[1],eqae[1],icodes[1],icaes[1]},type=numeric):
 Error, (in dsolve/numeric/DAE/checkconstraints) the initial conditions do not satisfy the algebraic constraints   error = .447e9, tolerance = .880e4, constraint = -2000000*(-1+y1(t))*exp(19.46229155000000000000*z1(t)-8.174162450000000000000)-2000000*y1(t)*exp(-19.46229155000000000000*z1(t)+8.174162450000000000000)+exp(38.92458310000000000000*z1(t)-11.79414868000000000000)-exp(-38.92458310000000000000*z1(t)+11.79414868000000000000)-100000
 > pars:=[mu=0.00001,q=1000,tint=1,tf=5001];
 (24)
 > Xexplicit:=2;
 (25)
 > ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));
 (26)
 > NODE:=nops(eqode):NAE:=nops(eqae);
 (27)
 > for XX from 1 to NODE do EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do;
 (28)
 > for XX from 1 to NAE do EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do;
 (29)
 > Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};
 (30)
 > Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};
 (31)
 > icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);
 (32)
 > for j from 1 to NAE do
 > EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);
 > end do;
 (33)
 > Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};
 (34)
 > if Xexplicit=1 then
 > sol:=dsolve(Sys,numeric,maxfun=0):
 > else
 > sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):
 > end if:
 >
 > for XX from 1 to NODE do a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red); end do:
 > for XX from NODE+1 to NODE+NAE do a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue); end do:
 > b1:=odeplot(sol,[t,y1(t)],0..1,color=red): b2:=odeplot(sol,[t,z1(t)],0..1,color=blue):
 > display(b1,b2,axes=boxed);
 > display(seq(a||x,x=1..NODE+NAE),axes=boxed);

End Example 2

 >

Example 3

 > restart;
 > with(plots):
 > eq1:=diff(y1(t),t)^2+diff(y1(t),t)*(y1(t)+1)+y1(t)=cos(diff(y1(t),t));
 (35)
 > solx:=dsolve({eq1,y1(0)=0},numeric):
 Error, (in dsolve/numeric/make_proc) Could not convert to an explicit first order system due to 'RootOf'
 > eqode:=[diff(y1(t),t)=z1(t)];
 (36)
 > eqae:=[subs(eqode,eq1)];
 (37)
 > icodes:=[y1(0)=0.0];
 (38)
 > icaes:=[z1(0)=0.0];
 (39)
 > pars:=[mu=0.1,q=1000,tint=1,tf=4];
 (40)
 > Xexplicit:=2;
 (41)
 > ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));
 (42)
 > NODE:=nops(eqode);NAE:=nops(eqae);
 (43)
 > for XX from 1 to NODE do EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do;
 (44)
 > for XX from 1 to NAE do EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do;
 (45)
 >
 > Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};
 (46)
 > Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};
 (47)
 > icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE);
 (48)
 > for j from 1 to NAE do
 > EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);
 > end do;
 (49)
 > Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};
 (50)
 > if Xexplicit=1 then
 > sol:=dsolve(Sys,numeric):
 > else
 > sol:=dsolve(Sys,numeric,stiff=true,implicit=true):
 > end if:
 >
 > for XX from 1 to NODE do a||XX:=odeplot(sol,[t,y||XX(t)],0..subs(pars,tf),color=red); end do:
 > for XX from NODE+1 to NODE+NAE do a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],0..subs(pars,tf),color=blue); end do:
 > display(seq(a||x,x=1..NODE+NAE),axes=boxed);

End Example 3

 >

Example 4

 > restart;
 > with(plots):
 > N:=11:h:=1/(N+1):
 > for i from 2 to N+1 do eq1[i]:=diff(y||i(t),t)=(y||(i+1)(t)-2*y||i(t)+y||(i-1)(t))/h^2-y||i(t)*(1+z||i(t));od:
 > for i from 2 to N+1 do eq2[i]:=0=(z||(i+1)(t)-2*z||i(t)+z||(i-1)(t))/h^2-(1-y||i(t)^2)*(exp(-z||i(t)));od:
 > eq1[1]:=(3*y1(t)-4*y2(t)+y3(t))/(2*h)=0: eq1[N+2]:=y||(N+2)(t)-1=0:
 > eq2[1]:=(3*z1(t)-4*z2(t)+1*z3(t))/(2*h)=0: eq2[N+2]:=z||(N+2)(t)=0:
 > eq1[1]:=subs(y1(t)=z||(N+3)(t),eq1[1]):
 > eq1[N+2]:=subs(y||(N+2)(t)=z||(N+4)(t),eq1[N+2]):
 > eqode:=[seq(subs(y1(t)=z||(N+3)(t),y||(N+2)(t)=z||(N+4)(t),eq1[i]),i=2..N+1)]:
 > eqae:=[eq1[1],eq1[N+2],seq(eq2[i],i=1..N+2)]:
 > icodes:=[seq(y||j(0)=1,j=2..N+1)]:
 > icaes:=[seq(z||j(0)=0,j=1..N+2),z||(N+3)(0)=1,z||(N+4)(0)=1]:
 > pars:=[mu=0.00001,q=1000,tint=1,tf=2]:
 > Xexplicit:=2:
 > ff:=subs(pars,1/2+1/2*tanh(q*(t-tint))):
 > NODE:=nops(eqode):NAE:=nops(eqae):
 > for XX from 1 to NODE do
 > EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:
 > for XX from 1 to NAE do
 > EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do:
 > Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)}:
 > Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)}:
 > icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE):
 > for j from 1 to NAE do
 > EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j):
 > end do:
 > Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:
 > if Xexplicit=1 then
 > sol:=dsolve(Sys,numeric,maxfun=0):
 > else
 > sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):
 > end if:
 >
 > for XX from 1 to NODE do
 > a||XX:=odeplot(sol,[t,y||(XX+1)(t)],1..subs(pars,tf),color=red): end do:
 > for XX from NODE+1 to NODE+NAE do
 > a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],1..subs(pars,tf),color=blue): end do:
 > display(seq(a||x,x=1..NODE),a||(NODE+NAE-1),a||(NODE+NAE),axes=boxed);

End of Example 4

 >

Sometimes the parameters of the switch function and perturbation need to be tuned to obtain propoer convergence. Below is Example 1 shown for several cases using the 'parameters' option in Maple's dsolve to compare how tuning parameters affects the solution

 > restart:
 > with(plots):
 > eqode:=[diff(y1(t),t)=-y1(t)^2+z1(t)]: eqae:=[cos(y1(t))-z1(t)^0.5=0]:
 > icodes:=[y1(0)=0.25]: icaes:=[z1(0)=0.8]:
 > pars:=[tf=5]:
 > Xexplicit:=2;
 (51)
 > ff:=subs(pars,1/2+1/2*tanh(q*(t-tint))):
 > NODE:=nops(eqode):NAE:=nops(eqae):
 > for XX from 1 to NODE do EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:
 > for XX from 1 to NAE do EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do:
 >
 > Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)}:
 > Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)}:
 > icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE):
 > for j from 1 to NAE do
 > EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j):
 > end do:
 > Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:
 > if Xexplicit=1 then
 > sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],maxfun=0):
 > else
 > sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],stiff=true,implicit=true):
 > end if:
 >
 > sol('parameters'=[0.1,1000,1]):
 > plot1:=odeplot(sol,[t-1,y1(t)],0..4,color=red): plot2:=odeplot(sol,[t-1,z1(t)],0..4,color=blue):
 > sol('parameters'=[0.001,10,1]):
 > plot3:=odeplot(sol,[t-1,y1(t)],0..4,color=red,linestyle=dot): plot4:=odeplot(sol,[t-1,z1(t)],0..4,color=blue,linestyle=dot):
 > display(plot1,plot2,plot3,plot4,axes=boxed);

In general, one has to decrease mu, and increase q and tint until convergence (example at t=3)

 > sol('parameters'=[0.001,10,1]):sol(3+1);
 (52)
 > sol('parameters'=[0.0001,100,10]):sol(3+10);
 (53)
 > sol('parameters'=[0.00001,1000,20]):sol(3+20);
 (54)
 >

The results have converged to 4 digits after the decimal. Of course, absolute and relative tolerances of the solvers can be modified if needed

## Combinations of numerical solutions from dsolve...

Hi! I am trying to plot and store in memory some specific combinations of the solutions of the systems of ODEs that I get numerically from dsolve for a particular range of the independent variable.

A particular case for my problem is the following system of stiff ODEs for two unknown functions f[0,0](x) and f[1,0](x) beween xini (where the Initial conditions are defined) and xfin, an arbitrary value of x. Note that rosebrock method does not work, and I can only solve it with lsode[adamsfull] or lsode[backfull]. I am attaching a maple file that shows what I have done.

############## System of ODEs that needs to be solved ####################################

 (1)

 (2)

###################################################################################

 (3)

 (4)

 (5)

The approach in that file works, however I have a question regarding the efficiency of my method, since I plan to extend the system to many more ODEs besides just 2 and also extend the range to a larger xfin. In this method, since I define the function to plot in terms of f01 and f02, wich are procedures, does this mean that for each x on the grid for the plot(ftoplot,x=xini..xfin) maple actually computes the solutions f00(x) and f01(x) and then forms the ftoplot combination and plots that specific point? If the default sampling of my interval is, say 1000 points, does it mean that the way I wrote it I will have 1000 invocations of the dsolve procedure, for each x in the sample? I am not sure, it seems to me that is the case. This would imply that instead of advancing the solution at each step maple starts over again from xini. How could I just avoid this behavior and instead have access to the values of ftoplot(x) in the range xini to xfin stored from one invocation of dsolve?

The ideal scenario for me would be to have f[0,0](x) and f[0,1](x) stored as an interpolated function between xini and xfin from the solutions of one invocation of dsolve prior to defining ftoplot. Can this be achieved in principle? How? Remember, i have to use method=lsode and range is not accepted.

## Most Optimal dsolve method for the following syste...

Hi!

I am trying to solve a large nxl system of coupled differential equations. Maple seems to have trouble even for small n's so I wanted to know if anyone has any suggestions. Take the case of the following system of ODEs for my unknown functions f[0,0](x) and f[1,0](x).

ODEs:= {diff(f[0, 0](x), x)+2.*f[0, 0](x)/x^5+.5000000000*f[0, 0](x)/x = -15.58845727*sin(.5773502693*x)/x^2+140.2961154*sin(.5773502693*x)/x^4-81.*cos(.5773502693*x)/x^3, diff(f[1, 0](x), x)+6.*f[1, 0](x)/x^5+1.500000000*f[1, 0](x)/x-1.*f[0, 0](x)/x = -15.58845727*sin(.5773502693*x)/x^2+25.98076212*sin(.5773502693*x)*(1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4)/x^2+140.2961154*sin(.5773502693*x)/x^4-233.8268591*sin(.5773502693*x)*(1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4)/x^4-81.*cos(.5773502693*x)/x^3+135.*cos(.5773502693*x)*(1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4)/x^3-20.78460970*sin(.5773502693*x)/x^6+6.000000004*cos(.5773502693*x)/x^5+62.35382908*sin(.5773502693*x)/x^8-36.00000002*cos(.5773502693*x)/x^7, f[0, 0](.1) = 1.503497680, f[1, 0](.1) = -.5011660086}

Following Preben Alsholm's suggestion from my previous thread I am using lsode[adamsfull], since no other method i have tried worked for this problem. I am currently using:

and it seems to work. I am wondering if there is a way to optimize this, as I will be extending my problem to n and l much larger than order unity numbers, therefore my system will contain about 10^4-10^5 equations. Solving this symple system of 2 equations takes a bit less than a second, but still it takes some time for the processor on my MBP. I am affraid it will be a nightmare for the full problem. Whats the most optimal dsolve option for this kind of problem? Any ideas?

I have also attempted dverk78, rkf45,rosenbrock, lsode(without the adamsfull option), and all failed for this particular system. Errors were:

1. For rkf45: Error, (in f00) cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up

2. For dverk78: Error, (in Soldverk78) cannot evaluate the solution past .1, step size < hmin, problem may be singular or error tolerance may be too small

3. For rosenbrock: Error, (in dsolve/numeric/SC/firststep) unable to evaluate the partial derivatives of f(x,y) for stiff solution

4. For lsode without [adamsfull]: Error, (in Sollsode) an excessive amount of work (greater than mxstep) was done

5. For default method with stiff=true and inplicit=true options: Error, (in dsolve/numeric/SC/firststep) unable to evaluate the partial derivatives of f(x,y) for stiff solution

## how to solve differential algebraic equations...

Hi all,

I try many ways to solve the below DAEs, for example stiff=true, method = rkf45_dae, method=rosenbrock..., but my DAEs can't be solved. In different ways, there are different errors, such as "Warning, cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up", ""

Would you please show me how to solve the below DAEs and plot their solutions?

## system of ODEs HDEDGP...

Dear guys!

I want to solve this system:

> c := 1: RC := 0.03: h := 1: m := 0.3:

> a := H(z)^2/(1+z)-h^2*(m*(1+z)^2+2*RC/(1+z));

> ode1 := diff(H(z),z) = (H(z)^4/((1+z)^2)+m*(1+z)*h^2*H(z)^2-2*m^2*h^4*(1+z)^4)/(2*H(z)*a) - (2*H(z)*h^2*RC)/(a*(1+z)^2) + (2*RC*h^2*(1+z)/(c*H(z)^2*a))*(H(z)^2/(1+z)^2-(H(z)^2/(1+z)-m*h^2*(1+z)^2)^2/(4*RC*h^2))^(3/2):

> ode2 := diff(M(z),z)=3*M(z)/(1+z)-2*M(z)*diff(H(z),z)/H(z):

> sys := {ode1, ode2}:

## Pushing dsolve to its limits

by: Maple 15

And so with this provocative title, "pushing dsolve to its limits" I want to share some difficulties I've been having in doing just that. I'm looking at a dynamic system of 3 ODEs. The system has a continuum of stationary points along a line. For each point on the line, there exist a stable (center) manifold, also a line, such that the point may be approached from both directions. However, simulating the converging trajectory has proven difficult.

I have simulated as...

## dsolve,numeric,stiff=true...

Hello all,

Is any anyone know how good is maple numeric stiff ode solver (stiff=true) ? What i mean is how many stiff ode's it can handle at a given time. I am having a hardtime solving a system of ode's (say 300 stiff ode's). Is there any other numeric stiff ode solver in maple that can handle large number of ode's?

I used islode advanced numeric solver for stiff ode's and it also has some limitations on the number of ode's and moreover it is slow.

## Numerical Stability...

Hi all,

I am having a very hard time with numerical stability. I am solving system of ode's (7-coupled ode's) using dsolve(stiff) and then using spline function for interpolation and finally solving system of pde's (2-coupled pde's) using pdsolve for one time step and solving all again for the next step. the solution is not stable and it requires very fine/small time step. Is there any procedure/method to improve the stability?

What is the stability criteria of dsolve...