Featured Posts

For the past thirty years, I have used several mathematical packages for problem solving and graphing. It all started with spreadsheet software that really helped speedup calculations compared to calculators. As many people do, once I had one tool I then started looking for another that would offer even more capabilities and features. I tested several of the very early math software but none really did all that I wanted until I came across Maple while I was working at SPAR Aerospace in Canada. For me, the rest is history. As long as I had a copy of Maple, it was all that I needed.

On occasions when I did not have a copy of this amazing software, I resorted to spreadsheets once more to complete fairly large and complex projects involving large databases and large numbers of calculations, especially when performing What-If scenarios. One distinct disadvantage of using a spreadsheet was the cryptic form of equation writing. I had to divide one long equation into several sections in different cells and then add them all up, which clearly is not good for documentation of the calculations. It is also very confusing for other engineers to know what that equation is or what it does. The development of the full engineering spreadsheet took months to complete, debug and verify. During this process, when I had errors, it was often very difficult to track exactly where the problem was, making the debugging process time consuming and sometimes very frustrating.

Having worked with Maple before, I remembered how easy it was to enter equations in a very familiar, readable math format. The real power of this software is that it allows you to write the equation(s) anyway you like and solve for any given parameter, unlike spreadsheets where you have to solve the problem first, by hand, for the parameter you want and then get the spreadsheet to calculate the value. I remember one time a few years ago when I wrote nine or ten simultaneous differential equations all in symbolic form and asked Maple to calculate certain parameters in a fully symbolic form. To my utmost disbelief, the answer came back within few minutes. With results in hand, I was able to quickly finish my research, and the results were published at PCIM Europe 2005 in “Distributed Gate ESR and its Effect on Shoot Through Performance at the Die Level”. I would never have gotten the results I needed if I was using a spreadsheet.

Even with much simpler systems of equations, finding solutions with a paper and pencil was never an easy task for me. It took a very long time, and even then there was no guarantee that I did not make copying errors, accidentally leave out a term, or make a calculation error. After I found the correct solution, I then had the problem of plotting the results, which I often needed in 3-D. Plotting allowed much deeper insights into the interdependency of all the parameters and made it easy for me to concentrate on the important ones without wasting any time. I was very happy when I could pass all these tasks onto Maple, which could do them much faster and more reliably then I ever could. Maple is a software that allows me to go beyond routine engineering calculations and gives me the tools to reach levels of insight and understanding that were completely out of reach of the average engineer until a few years ago.

For the record, I have no business affiliations with Maplesoft. I’m writing this article because Maple makes such a difference in my work that I feel it is important to share my experiences so other engineers can get the same benefits.

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.

http://depts.washington.edu/maple/pubs/U_Apprach_FULL_DRAFT.pdf  

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)];

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];

eqae := [cos(y1(t))-z1(t)^.5 = 0]

(2)

Enter all initial conditions for differential variables in icodes

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

icodes := [y1(0) = .25]

(3)

Enter all intial conditions for algebraic variables in icaes

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

icaes := [z1(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];

pars := [mu = .1, q = 1000, tint = 1, tf = 5]

(5)

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

Xexplicit:=2;

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

ff := 1/2+(1/2)*tanh(1000*t-1000)

(7)

NODE:=nops(eqode);NAE:=nops(eqae);

NODE := 1

NAE := 1

(8)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff;
end do;

EQODE1 := diff(y1(t), t) = (-y1(t)^2+z1(t))*(1/2+(1/2)*tanh(1000*t-1000))

(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;

EQAE1 := -.1*sin(y1(t))*(diff(y1(t), t))-0.5e-1*(diff(z1(t), t))/z1(t)^.5 = -cos(y1(t))+z1(t)^.5

(10)

 

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(11)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(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);

icsn := y1(t) = .25, z1(t) = .8

(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)};

Sys := {-0.1824604200e-1-0.5590169945e-1*(diff(z1(t), t)) = -cos(y1(t))+z1(t)^.5, y1(0) = .25, z1(0) = .8, diff(y1(t), t) = (-y1(t)^2+z1(t))*(1/2+(1/2)*tanh(1000*t-1000))}

(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;

eq1 := diff(y1(t), t) = j1*W/(F*rho*V)

(15)

eq2:=j1+j2=iapp;

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

j1 := io1*(2*(1-y1(t))*exp(.5*F*(z1(t)-phi1)/(R*T))-2*y1(t)*exp(-.5*F*(z1(t)-phi1)/(R*T)))

(17)

j2:=io2*(exp((F/R/T)*(z1(t)-phi2))-exp((-F/R/T)*(z1(t)-phi2)));

j2 := io2*(exp(F*(z1(t)-phi2)/(R*T))-exp(-F*(z1(t)-phi2)/(R*T)))

(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};

params := {F = 96487, R = 8.314, T = 298.15, V = 0.1e-4, W = 92.7, io1 = 0.1e-3, io2 = 0.1e-9, rho = 3.4, iapp = 0.1e-4, phi1 = .420, phi2 = .303}

(19)

eqode:=[subs(params,eq1)];

eqode := [diff(y1(t), t) = 0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450)]

(20)

eqae:=[subs(params,eq2)];

eqae := [0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)+0.1e-9*exp(38.92458310*z1(t)-11.79414868)-0.1e-9*exp(-38.92458310*z1(t)+11.79414868) = 0.1e-4]

(21)

icodes:=[y1(0)=0.05];

icodes := [y1(0) = 0.5e-1]

(22)

icaes:=[z1(0)=0.7];

icaes := [z1(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];

pars := [mu = 0.1e-4, q = 1000, tint = 1, tf = 5001]

(24)

Xexplicit:=2;

Xexplicit := 2

(25)

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));

ff := 1/2+(1/2)*tanh(1000*t-1000)

(26)

NODE:=nops(eqode):NAE:=nops(eqae);

NAE := 1

(27)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff:
end do;

EQODE1 := diff(y1(t), t) = (0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450))*(1/2+(1/2)*tanh(1000*t-1000))

(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;

EQAE1 := -0.2e-8*(diff(y1(t), t))*exp(19.46229155*z1(t)-8.174162450)+0.3892458310e-7*(1-y1(t))*(diff(z1(t), t))*exp(19.46229155*z1(t)-8.174162450)-0.2e-8*(diff(y1(t), t))*exp(-19.46229155*z1(t)+8.174162450)+0.3892458310e-7*y1(t)*(diff(z1(t), t))*exp(-19.46229155*z1(t)+8.174162450)+0.3892458310e-13*(diff(z1(t), t))*exp(38.92458310*z1(t)-11.79414868)+0.3892458310e-13*(diff(z1(t), t))*exp(-38.92458310*z1(t)+11.79414868) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868)

(29)

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(30)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(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);

icsn := y1(t) = 0.5e-1, z1(t) = .7

(32)

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);

end do;

EQAEX1 := -0.2e-8*(0.5368903705e-2*exp(5.449441630)-0.2825738792e-3*exp(-5.449441630))*exp(5.449441630)+0.3697835394e-7*(diff(z1(t), t))*exp(5.449441630)-0.2e-8*(0.5368903705e-2*exp(5.449441630)-0.2825738792e-3*exp(-5.449441630))*exp(-5.449441630)+0.1946229155e-8*(diff(z1(t), t))*exp(-5.449441630)+0.3892458310e-13*(diff(z1(t), t))*exp(15.45305949)+0.3892458310e-13*(diff(z1(t), t))*exp(-15.45305949) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868)

(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)};

Sys := {-0.5810962488e-6+0.8802389238e-5*(diff(z1(t), t)) = 0.1e-4-0.2e-3*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)+0.2e-3*y1(t)*exp(-19.46229155*z1(t)+8.174162450)-0.1e-9*exp(38.92458310*z1(t)-11.79414868)+0.1e-9*exp(-38.92458310*z1(t)+11.79414868), y1(0) = 0.5e-1, z1(0) = .7, diff(y1(t), t) = (0.5651477584e-2*(1-y1(t))*exp(19.46229155*z1(t)-8.174162450)-0.5651477584e-2*y1(t)*exp(-19.46229155*z1(t)+8.174162450))*(1/2+(1/2)*tanh(1000*t-1000))}

(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));

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)];

eqode := [diff(y1(t), t) = z1(t)]

(36)

eqae:=[subs(eqode,eq1)];

eqae := [z1(t)^2+z1(t)*(y1(t)+1)+y1(t) = cos(z1(t))]

(37)

icodes:=[y1(0)=0.0];

icodes := [y1(0) = 0.]

(38)

icaes:=[z1(0)=0.0];

icaes := [z1(0) = 0.]

(39)

pars:=[mu=0.1,q=1000,tint=1,tf=4];

pars := [mu = .1, q = 1000, tint = 1, tf = 4]

(40)

Xexplicit:=2;

Xexplicit := 2

(41)

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint)));

ff := 1/2+(1/2)*tanh(1000*t-1000)

(42)

NODE:=nops(eqode);NAE:=nops(eqae);

NODE := 1

NAE := 1

(43)

for XX from 1 to NODE do
EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff:
end do;

EQODE1 := diff(y1(t), t) = z1(t)*(1/2+(1/2)*tanh(1000*t-1000))

(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;

EQAE1 := .1*sin(z1(t))*(diff(z1(t), t))+.2*z1(t)*(diff(z1(t), t))+.1*(diff(z1(t), t))*(y1(t)+1)+.1*z1(t)*(diff(y1(t), t))+.1*(diff(y1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t)

(45)

 

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)};

Dvars1 := {diff(z1(t), t) = D1}

(46)

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)};

Dvars2 := {D1 = diff(z1(t), t)}

(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);

icsn := y1(t) = 0., z1(t) = 0.

(48)

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j);

end do;

EQAEX1 := .1*sin(0.)*(diff(z1(t), t))+.1*(diff(z1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t)

(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)};

Sys := {.1*(diff(z1(t), t)) = cos(z1(t))-z1(t)^2-z1(t)*(y1(t)+1)-y1(t), y1(0) = 0., z1(0) = 0., diff(y1(t), t) = z1(t)*(1/2+(1/2)*tanh(1000*t-1000))}

(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;

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

[t = 4., y1(t) = .738587929442734, z1(t) = .546472878850096]

(52)

sol('parameters'=[0.0001,100,10]):sol(3+10);

[t = 13., y1(t) = .738684397167344, z1(t) = .546618936273638]

(53)

sol('parameters'=[0.00001,1000,20]):sol(3+20);

[t = 23., y1(t) = .738694113087217, z1(t) = .546633473784526]

(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

 

Download SolvingDAEsinMaple.mws


MaplePrimes Questions Recent Unanswered Maple MapleSim Maple T.A.