Appendix A
Maple code for Examples 14 from "Extending Explicit and Linealry Implicit ODE Solvers for Index1 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
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

(3) 
Enter all intial conditions for algebraic variables in icaes

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

(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 = .745e1, tolerance = .559e6, constraint = cos(y1(t))z1(t)^.5000000000000000000000


> 
ff:=subs(pars,1/2+1/2*tanh(q*(ttint)));


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


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


(9) 
> 
for XX from 1 to NAE do EQAEXX:=subs(pars,mu*(diff(rhs(eqae[XX])lhs(eqae[XX]),t))=rhs(eqae[XX])lhs(eqae[XX])); end do;


(10) 
> 
Dvars1:={seq(diff(zx(t),t)=Dx,x=1..NAE)};


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


(12) 
> 
icsn:=seq(subs(yx(0)=yx(t),icodes[x]),x=1..NODE),seq(subs(zx(0)=zx(t),icaes[x]),x=1..NAE);


(13) 
> 
EQAEXj:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAEj))=rhs(EQAEj);

> 
Sys:={seq(EQODEx,x=1..NODE),seq(EQAEXx,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};


(14) 
> 
sol:=dsolve(Sys,numeric):

> 
sol:=dsolve(Sys,numeric,stiff=true,implicit=true): end if:

> 
for XX from 1 to NODE do aXX:=odeplot(sol,[t,yXX(t)],0..subs(pars,tf),color=red); end do:

> 
for XX from NODE+1 to NODE+NAE do aXX:=odeplot(sol,[t,z(XXNODE)(t)],0..subs(pars,tf),color=blue); end do:

> 
display(seq(ax,x=1..NODE+NAE),axes=boxed);

End Example 1
Example 2
> 
eq1:=diff(y1(t),t)=j1*W/F/rho/V;


(15) 

(16) 
> 
j1:=io1*(2*(1y1(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=1e5,io1=1e4,io2=1e10,iapp=1e5,rho=3.4};


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


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


(21) 

(22) 

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

(25) 
> 
ff:=subs(pars,1/2+1/2*tanh(q*(ttint)));


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


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


(28) 
> 
for XX from 1 to NAE do EQAEXX:=subs(pars,mu*(diff(rhs(eqae[XX])lhs(eqae[XX]),t))=rhs(eqae[XX])lhs(eqae[XX])); end do;


(29) 
> 
Dvars1:={seq(diff(zx(t),t)=Dx,x=1..NAE)};


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


(31) 
> 
icsn:=seq(subs(yx(0)=yx(t),icodes[x]),x=1..NODE),seq(subs(zx(0)=zx(t),icaes[x]),x=1..NAE);


(32) 
> 
EQAEXj:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAEj))=rhs(EQAEj);


(33) 
> 
Sys:={seq(EQODEx,x=1..NODE),seq(EQAEXx,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};


(34) 
> 
sol:=dsolve(Sys,numeric,maxfun=0):

> 
sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):

> 
for XX from 1 to NODE do aXX:=odeplot(sol,[t,yXX(t)],0..subs(pars,tf),color=red); end do:

> 
for XX from NODE+1 to NODE+NAE do aXX:=odeplot(sol,[t,z(XXNODE)(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(ax,x=1..NODE+NAE),axes=boxed);

End Example 2
Example 3
> 
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) 

(38) 

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


(40) 

(41) 
> 
ff:=subs(pars,1/2+1/2*tanh(q*(ttint)));


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


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


(44) 
> 
for XX from 1 to NAE do EQAEXX:=subs(pars,mu*(diff(rhs(eqae[XX])lhs(eqae[XX]),t))=rhs(eqae[XX])lhs(eqae[XX])); end do;


(45) 
> 
Dvars1:={seq(diff(zx(t),t)=Dx,x=1..NAE)};


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


(47) 
> 
icsn:=seq(subs(yx(0)=yx(t),icodes[x]),x=1..NODE),seq(subs(zx(0)=zx(t),icaes[x]),x=1..NAE);


(48) 
> 
EQAEXj:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAEj))=rhs(EQAEj);


(49) 
> 
Sys:={seq(EQODEx,x=1..NODE),seq(EQAEXx,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)};


(50) 
> 
sol:=dsolve(Sys,numeric):

> 
sol:=dsolve(Sys,numeric,stiff=true,implicit=true):

> 
for XX from 1 to NODE do aXX:=odeplot(sol,[t,yXX(t)],0..subs(pars,tf),color=red); end do:

> 
for XX from NODE+1 to NODE+NAE do aXX:=odeplot(sol,[t,z(XXNODE)(t)],0..subs(pars,tf),color=blue); end do:

> 
display(seq(ax,x=1..NODE+NAE),axes=boxed);

End Example 3
Example 4
> 
for i from 2 to N+1 do eq1[i]:=diff(yi(t),t)=(y(i+1)(t)2*yi(t)+y(i1)(t))/h^2yi(t)*(1+zi(t));od:

> 
for i from 2 to N+1 do eq2[i]:=0=(z(i+1)(t)2*zi(t)+z(i1)(t))/h^2(1yi(t)^2)*(exp(zi(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(yj(0)=1,j=2..N+1)]:

> 
icaes:=[seq(zj(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]:

> 
ff:=subs(pars,1/2+1/2*tanh(q*(ttint))):

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

> 
for XX from 1 to NODE do

> 
EQODEXX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:

> 
for XX from 1 to NAE do

> 
EQAEXX:=subs(pars,mu*(diff(rhs(eqae[XX])lhs(eqae[XX]),t))=rhs(eqae[XX])lhs(eqae[XX])); end do:

> 
Dvars1:={seq(diff(zx(t),t)=Dx,x=1..NAE)}:

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

> 
icsn:=seq(subs(yx(0)=yx(t),icodes[x]),x=1..NODE),seq(subs(zx(0)=zx(t),icaes[x]),x=1..NAE):

> 
EQAEXj:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAEj))=rhs(EQAEj):

> 
Sys:={seq(EQODEx,x=1..NODE),seq(EQAEXx,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:

> 
sol:=dsolve(Sys,numeric,maxfun=0):

> 
sol:=dsolve(Sys,numeric,stiff=true,implicit=true,maxfun=0):

> 
for XX from 1 to NODE do

> 
aXX:=odeplot(sol,[t,y(XX+1)(t)],1..subs(pars,tf),color=red): end do:

> 
for XX from NODE+1 to NODE+NAE do

> 
aXX:=odeplot(sol,[t,z(XXNODE)(t)],1..subs(pars,tf),color=blue): end do:

> 
display(seq(ax,x=1..NODE),a(NODE+NAE1),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
> 
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]:


(51) 
> 
ff:=subs(pars,1/2+1/2*tanh(q*(ttint))):

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

> 
for XX from 1 to NODE do EQODEXX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:

> 
for XX from 1 to NAE do EQAEXX:=subs(pars,mu*(diff(rhs(eqae[XX])lhs(eqae[XX]),t))=rhs(eqae[XX])lhs(eqae[XX])); end do:

> 
Dvars1:={seq(diff(zx(t),t)=Dx,x=1..NAE)}:

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

> 
icsn:=seq(subs(yx(0)=yx(t),icodes[x]),x=1..NODE),seq(subs(zx(0)=zx(t),icaes[x]),x=1..NAE):

> 
EQAEXj:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAEj))=rhs(EQAEj):

> 
Sys:={seq(EQODEx,x=1..NODE),seq(EQAEXx,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:

> 
sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],maxfun=0):

> 
sol:=dsolve(Sys,numeric,'parameters'=[mu,q,tint],stiff=true,implicit=true):

> 
sol('parameters'=[0.1,1000,1]):

> 
plot1:=odeplot(sol,[t1,y1(t)],0..4,color=red): plot2:=odeplot(sol,[t1,z1(t)],0..4,color=blue):

> 
sol('parameters'=[0.001,10,1]):

> 
plot3:=odeplot(sol,[t1,y1(t)],0..4,color=red,linestyle=dot): plot4:=odeplot(sol,[t1,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
