vsubramanian - Posts - MaplePrimes

Maple helps in simulating multiscale-mul...

Posted: Dr. Venkat Subramanian 431 Product: Maple

October 06 2020

1 0

We recently published a paper on multiscale-multidomain simulation of battery models.

https://iopscience.iop.org/article/10.1149/1945-7111/abb37b

Some challenges are listed at

https://twitter.com/UT_MAPLE/status/1311356957941522438

Maple and symbolic math can play a critical role in solving many challenging problems. For example, consider a seemingly-trivial problem

uxx+uyy = 0

x = 0 and x = 1, ux = 0 for all y

y = 1, u = 0, for all x

y = 0, u =1, 0<=x<=0.5

y = 0, uy = 0, 0.5<x<=1

There is a singularity at (0.5,0) and most numerical methods will have trouble there. In these equations, uxx means the second derivative of u with respect to x.

Maple can help solve this problem with conformal mapping to achieve arbitrary precision. As of today, machine precision is not possible with any numerical method even with millions of Degrees of Freedom. The Maple code is given below. A FEM code is given below as well. Models like this can benefit from Maple adding parallel sparse direct and iterative solvers.

>

restart;

>	Digits:=15;

(1)

The domain is tranformed from Z = X+IY to w. The points tranformed are

>	Zdomain:=[[0,0],[1,0],[1,1],[0,1]];

(2)

>	wdomain:=[0,1,1+a,a+2];

(3)

>	eq1:=diff(Z(w),w)=-I*K1/sqrt(w)/sqrt(w-1)/sqrt(w-a-1)/sqrt(w-a-2);

(4)

a is not known apriori. The value of a should be found to make sure [1,1] in the Z coordinate is transformed to [1+a] in the w coordinate.

>	a:=sqrt(2)-1;

(5)

Value of K1 is found using the transformation of [1,0] to 1 in the w coordinate

>	eq11:=1=int(rhs(eq1),w=0..1);

(6)

>	K1:=solve(eq11,K1);

K1 := -(1/2)*(2^(1/2)/EllipticK(2^(1/2)/2))

(7)

The height Y in the Z coordinate is found by integrating from 1 to 1+a in the w coordinate.

>	simplify(int(rhs(eq1),w=1..1+a));

(2*I)*2^(1/2)*EllipticK((2^(1/2)-1)/(2^(1/2)+1))/(EllipticK(2^(1/2)/2)*(2^(1/2)+1))

(8)

>

evalf(%);

(9)

The choice of a = sqrt(2)-1 gives the height of 1 for Z coordinate.

>	eval(simplify(int(rhs(eq1),w=0..1+a)));

(2^(1/2)*EllipticK(2^(1/2)/2)+EllipticK(2^(1/2)/2)+(2*I)*2^(1/2)*EllipticK((2^(1/2)-1)/(2^(1/2)+1)))/(EllipticK(2^(1/2)/2)*(2^(1/2)+1))

(10)

>

evalf(%);

(11)

integration from 0 to 1+a in the w coordinate gives 1,1 in the Z coordinate.

>	simplify(int(rhs(eq1),w=0..1/sqrt(2)));

EllipticF(2^(1/2)/(2+2^(1/2))^(1/2), 2^(1/2)/2)/EllipticK(2^(1/2)/2)

(12)

>

evalf(%);

(13)

Integrating w from 0 to wmid =1/sqrt(2) gives the point 0.5,0 in the Z coordinate

>	wmid:=1/sqrt(2);

(14)

Next w domain is transformed to Z2 domain Z2 = X2+IY2

>	Z2domain:=[[0,0],[1,0],[1,H],[0,H]];

(15)

>	wdomain2:=[0,1/sqrt(2),1+a,a+2];

(16)

>	eq2:=diff(Z2(w),w)=-I*K2/sqrt(w)/sqrt(w-wmid)/sqrt(w-1-a)/sqrt(w-a-2);

eq2 := diff(Z2(w), w) = -(2*I)*K2/(w^(1/2)*(4*w-2*2^(1/2))^(1/2)*(w-2^(1/2))^(1/2)*(w-2^(1/2)-1)^(1/2))

(17)

K2 is found based on the transformation of wmid to [1,0] in Z2 coordinate.

>	eq21:=1=int(rhs(eq2),w=0..wmid);

eq21 := 1 = -2*K2*EllipticK(1/(((2+2^(1/2))^(1/2))))/(2^(1/2)+1)^(1/2)

(18)

>	K2:=solve(eq21,K2);

K2 := -(1/2)*((2^(1/2)+1)^(1/2)/EllipticK(1/(((2+2^(1/2))^(1/2)))))

(19)

The corner 0,1 in the Z coordinate is mapped by integrating eq2 from 0 to 1 in the w coordinate

>	int(rhs(eq2),w=0..1);

1+EllipticF((2-2^(1/2))^(1/2), ((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*I/EllipticK(1/(((2+2^(1/2))^(1/2))))

(20)

>	corner:=evalf(%);

(21)

This is the point 1,0 in the original coordinate.

The height in the Z2 coorinate is found by integrating eq2 from wmid to 1+a.

>	ytot:=int(rhs(eq2),w=wmid..1+a);

ytot := EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*I/EllipticK(1/(((2+2^(1/2))^(1/2))))

(22)

>

evalf(%);

(23)

The magnitude in the Y direction is given by the coefficient of I, the imaginary number

>	Ytot:=coeff(ytot,I);

Ytot := EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))/EllipticK(1/(((2+2^(1/2))^(1/2))))

(24)

The analytial solution in the Z2 corordinate is a line in Y2 to satisfy simple zero flux conditions at X2 = 0 and X2 = 1.

>	phianal:=1+b*Y2;

(25)

The value of phi is zero at Y2 = Ytot (originally the cathode domain in the Z domain)

>	bc:=subs(Y2=Ytot,phianal)=0;

bc := 1+b*EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))/EllipticK(1/(((2+2^(1/2))^(1/2)))) = 0

(26)

>	b:=solve(bc,b);

b := -EllipticK(1/(((2+2^(1/2))^(1/2))))/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(27)

The analytical solution is given by

>

phianal;

1-EllipticK(1/(((2+2^(1/2))^(1/2))))*Y2/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(28)

>	evalf(phianal);

(29)

The potential at the corner is given by substituting the imaginary value of corner for Y2 in phinanal)

>	phicorner:=subs(Y2=Im(corner),phianal);

phicorner := 1-.559419351518322*EllipticK(1/(((2+2^(1/2))^(1/2))))/EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))

(30)

>	evalf(phicorner);

(31)

local flux/current density calculation, written in terms of w is

>	curr:=b*rhs(eq2)/rhs(eq1);

curr := -(2^(1/2)+1)^(1/2)*2^(1/2)*EllipticK(2^(1/2)/2)*(w-1)^(1/2)/(EllipticK(((2^(1/2)+1)/(2+2^(1/2)))^(1/2))*(4*w-2*2^(1/2))^(1/2))

(32)

average flux/current density calculation for the anode

>	currave:=int((curr),w=0..wmid)/wmid;

currave := -(1/2)*((2^(1/2)+1)^(1/2)*(2^(1/2)-1)*EllipticK(2^(1/2)/2)*(2^(3/4)+2^(1/4)+arctanh(2^(3/4)/2))*2^(1/2)/EllipticK(2^(3/4)/2))

(33)

>	Digits:=25:

The average current density at Y =0, local current density at X = 0,Y=0 and potential at X=1,Y=0 (Corner) can be used to study convergence of FEM and other numerical methods

>	evalf(currave),evalf(subs(w=0,curr)),evalf(phicorner);

(34)

>

Download Conformalmapping.mws

This FEM code is for solving Laplace's equation with primary current distribution considered in Model 1.
This code is based on FEM weak-form. Biquadratic Lagrange shape functions (9nodes in an element) are used.

>

restart;

>	with(LinearAlgebra):

>	Lx:=1: #length in X

>	Ly:=1: #length in Y

>	nx:=10: #number of elements in X (even numbers only)

>	ny:=10: #number of elements in Y, to be kept same as nx in this version

>	hx:=Lx/nx: #element size x

>	hy:=Ly/ny: #element size y

Procedure to perform numerical integration on shape functions to obtain local matrices (can be replaced with analytical integration for this particular problem)
-Shape functions are also used as weight functions in applying weak formulation. Numerical integration is done using Simpson's rule.
-Local cartesian coordinates x,y are converted to natural coordinates zeta and eta. This transformation is not required for this simple geomerty but useful in general. zeta and eta are obtained by scaling x and y with hx/2 and hy/2, respectively, in this code.

>	localmatrices:=proc(a1,a2,a3,b1,b2,b3,q1,q2) global Kx,Ky,Nx,Ny,zeta,eta,c; local A,dAdzeta,dAdeta,y,x,J,terms,i,j,k,l,dx,dy,fx,fy,fxy,fyy,dzeta,deta,J1,J2;

>	A:=[(1-zeta)zeta(1-eta)eta/4,-(1-(zeta)^2)(1-eta)eta/2,-(1+zeta)zeta(1-eta)eta/4,(1-(eta)^2)(1+zeta)zeta/2,(1+zeta)zeta(1+eta)eta/4,(1-(zeta)^2)(1+eta)eta/2,-(1-zeta)zeta(1+eta)eta/4,-(1-(eta)^2)(1-zeta)(zeta)/2,(1-(zeta)^2)*(1-(eta)^2)]; #bi quadratic langrange shape functions

>	dAdzeta:=diff(A,zeta);

>	dAdeta:=diff(A,eta);

>	y:=[-a1,-a2,-a3,0,a3,a2,a1,0,0];x:=[-b1,0,b1,b2,b3,0,-b3,-b2,0];

>	J:=simplify(add(dAdzeta[i]x[i],i=1..9)add(dAdeta[i]y[i],i=1..9)-add(dAdzeta[i]y[i],i=1..9)add(dAdeta[i]x[i],i=1..9));

>	Nx:=[seq((simplify(add(dAdeta[i]y[i],i=1..9))dAdzeta[j]-simplify(add(dAdzeta[i]y[i],i=1..9))dAdeta[j])/simplify(J),j=1..9)];

>	Ny:=[seq((-dAdzeta[j]simplify(add(dAdeta[i]x[i],i=1..9))+dAdeta[j]simplify(add(dAdzeta[i]x[i],i=1..9)))/simplify(J),j=1..9)];

>	Kx:=Matrix(nops(A),nops(A),datatype=float[8]):

>	Ky:=Matrix(nops(A),nops(A),datatype=float[8]): c:=Vector(nops(A),datatype=float[8]):

>	terms:=20:#number of terms for numerical integration dzeta:=2/terms: deta:=2/terms:

>

for i from 1 to nops(Nx) do #loop to obtain local matrices

for j from 1 to nops(Ny) do
Kx[i,j]:=0;
Ky[i,j]:=0;

for k from 0 to terms do #outer loop double integration, integration in zeta

if k = 0 then fx[k]:= subs(zeta=-1,Nx[i]*Nx[j]*J); fy[k]:= subs(zeta=-1,Ny[i]*Ny[j]*J);
elif k = terms then
fx[k]:= subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:= subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J);
elif irem(k,2) = 0 then
fx[k]:= 2*subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:=     2*subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J);
else fx[k]:= 4*subs(zeta=-1+(k*dzeta),Nx[i]*Nx[j]*J);
fy[k]:=     4*subs(zeta=-1+(k*dzeta),Ny[i]*Ny[j]*J); end if;

for l from 0 to terms do #inner loop double integration, integration in eta

if l = 0 then fxy[l]:= subs(eta=-1,fx[k]); fyy[l]:= subs(eta=-1,fy[k]);
elif l = terms then fxy[l]:= subs(eta=-1+(l*deta),fx[k]); fyy[l]:= subs(eta=-1+(l*deta),fy[k]);
elif irem(l,2) = 0 then fxy[l]:= 2*subs(eta=-1+(l*deta),fx[k]); fyy[l]:= 2*subs(eta=-1+(l*deta),fy[k]);
else fxy[l]:=4*subs(eta=-1+(l*deta),fx[k]); fyy[l]:=4*subs(eta=-1+(l*deta),fy[k]); end if;
Kx[i,j]:=Kx[i,j]+fxy[l];
Ky[i,j]:=Ky[i,j]+fyy[l];

end do;

end do;
Kx[i,j]:=Kx[i,j]*dzeta*deta/9;
Ky[i,j]:=Ky[i,j]*dzeta*deta/9;
end do;
end do:

>

end proc:

>	n:=nx*ny; #total number of elements

(1)

>	Nx1:=nx*2+1: #number of nodes in x in one row

>	N:=Nx1(2ny+1); # total number of nodes/equations

(2)

>	K:=Matrix(N,N,storage=sparse): # global K matrix

>	C:=Vector(N,storage=sparse): # global c matrix

>	L2G:=Matrix(n,9): #mapping matrix - each row has node numbers for each element

>	l:=1:k:=1: localmatrices(hy/2,hy/2,hy/2,hx/2,hx/2,hx/2,0,0): kx:=copy(Kx):ky:=copy(Ky):c0:=copy(c):

>

for i from 1 to n do #modifying,adding and assembling matrices to get global matrix

if i<=nx/2 then
a1:=copy(kx); a2:=copy(ky); a3:=0; a1[1..3,1..9]:=IdentityMatrix(3,9); a2[1..3,1..9]:=Matrix(3,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[1..3]:=1.0;

elif i=nx/2+1 then a1:=copy(kx); a2:=copy(ky); a3:=0; a1[1,1..9]:=IdentityMatrix(1,9); a2[1,1..9]:=Matrix(1,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[1]:=1.0;

elif i>nx*(ny-1) then a1:=copy(kx); a2:=copy(ky); a3:=0; a1[5..7,5..7]:=IdentityMatrix(3,3);
a1[5..7,1..4]:=ZeroMatrix(3,4);
a1[5..7,8..9]:=ZeroMatrix(3,2); a2[5..7,1..9]:=Matrix(3,9,shape=zero); a4:=a1+a2; c:=copy(c0); c[5..7]:=0;

else a1:=kx; a2:=ky; a3:=0;a4:=a1+a2; c:=c0; end if;

L2G[i,1..9]:=Matrix([l,l+1,l+2,l+2+Nx1,l+2+Nx1*2,l+1+Nx1*2,l+Nx1*2,l+Nx1,l+1+Nx1]):

k:=k+1:

if k>nx then k:=1; l:=l+Nx1+3;

else l:=l+2; end if:

indx2:=L2G[i,1..9]:
indx2:=convert(indx2,list):

C[indx2]:=C[indx2]+c[1..9];
c[1..9];

for i1 from 1 to 9 do
indx1:=L2G[i,i1]:
K[indx1,indx2]:=K[indx1,indx2]+a4[i1,1..9]:
end do:

end do:

>	phi:=LinearSolve(K,C,method=SparseLU): #linear set of equations solved using Sparse LU solver

>	phi_at(1,0):=phi[Nx1];

(3)

>	dNdy:=copy(Ny):

>	dNdy_bottom_left:=subs(eta=-1,zeta=-1,dNdy):

>	current_at(0,0):=add(dNdy_bottom_left[i]*phi[L2G[1,i]],i=1..nops(Ny));

(4)

>	if irem(nx/2,2)=0 then current_at(0,0.25):=add(dNdy_bottom_left[i]phi[L2G[nx/4+1,i]],i=1..nops(Ny)); else dNdy_bottom_center:=subs(eta=-1,zeta=0,dNdy): current_at(0,0.25):=add(dNdy_bottom_center[i]phi[L2G[(nx/2+1)/2,i]],i=1..nops(Ny)); end if;

(5)

>

Download FEM_2D.mws

A new and efficient code in Maple for Bo...

Posted: Dr. Venkat Subramanian 431 Product: Maple

August 23 2017

5 10

A new code based on higher derivative method has been implemented in Maple. A sample code is given below and explained. Because of the symbolic nature of Maple, this method works very well for a wide range of BVP problems.

The code solves BVPs written in the first order form dy/dx = f (Maple’s dsolve numeric converts general BVPs to this form and solves).

The code can handle unknown parameters in the model if sufficient boundary conditions are provided.

This code has been tested from Maple 8 to Maple 2017. For Digits:=15 or less, this code works in all of the Maple versions tested.

Most problems can be solved with Digits:=15 with atol = 1e-10 or so. This code can be used to get a tolerance value of 1e-20 or any high precision as needed by changing the number of Digits accordingly. This may be needed if the original variables are not properly sacled. With arbitrarily high Digits, the code fails in Maple 18 or later version, etc because Maple does not support SparseDirect Solver at high precision in some of the versions (hopefuly this bug can be removed in the future versions).

For simple problems, Maple’s dsolve/numeric is superior to the code developed as it is implemented in hardware floats. For large scale problems and stiff problems, the method developed is much more superior to Maple and comparable to (and often times better than) state of the art codes for BVPs - bvp4c (MATLAB), COLSYS,TWPBVP, etc.

The code, as written, cannot be used for problems with a singularity at end points (doable in the future). In addition, mixed boundary conditions are not supported in this version of the code (for example, y1(1)=y2(0)). Future updates will include the application of this approach for DAE-BVPs, currently not supported by Maple’s dsolve/numeric command.

A paper has been submitted to JCAM. I welcome feedback on the code and solicit input from Mapleprimes members if they are able to test (and break this code) for any BVP.

PDF of the paper submitted, example maple code and the solver as a text file needed are uploaded here. Additional examples are hosted on my website at http://depts.washington.edu/maple/HDM.html

>

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

Troesch's problem
This is an inherently unstable, difficult, nonlinear, two-point BVP formulated by Weibel and Troesch that describes the confinement of a place column by radiation pressure. Increasing epsilon increases the stiffness of the ODE.
1. E.S. Weibel, On the confinement of a plasma by magnetostatic fields, Phys. Fluids. 2 (1959) 52-56.
2. B. Troesch, A simple approach to a sensitive two-point boundary value problem, J. Comput. Phys. 21 (1976) 279-290.

Introduction
The package HDM solves boundary value problems (BVPs) using higher derivative methods (HDM) in Maple®. We explain how to solve BVPs using this package. HDM can numerically solve BVPs of ordinary differential equations (ODEs) of the form shown is the fowllowing example.

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

>

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;

(1)

>

Enter the first-order ODEs into EqODEs list.

>	EqODEs:=[diff(y1(x),x)=y2(x),diff(y2(x),x)=epsilonsinh(epsilony1(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)]);

(3)

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

(4)

>

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

>	Range:=[0.,1.];

(5)

>

List any known parameters in the list.

>	pars:=[epsilon=2];

(6)

>

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

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

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

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

>

Present some details of the solution.

>	sol[4]; # final derivative

(10)

>	sol[5]; # Maximum local RMSE

(11)

>

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

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

(12)

>

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

(13)

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

>

Download Example_3_Troesch.mws

Please bring back dsolve numeric method ...

Posted: Dr. Venkat Subramanian 431 Product: Maple

July 20 2016

6 5

In a recent conversation I explained whyLSODE was giving wrong results (http://www.mapleprimes.com/questions/210948-Can-We-Trust-Maple#comment230167). After a lot of confusions and weird infinite loops for answers, it turned out that Newton Raphson was not properly done.

Both LSODE and MEBDFI are currently incompletely implemented (only one iteration is done instead of Newton Raphson till convergence). Maplesoft should update the help files accordingly.

The post below explains how better results are obtained with method = mgear. To run the command mgear you will need Maple 6 or earlier versions. For lsode, any current version is fine. Unfortunately Maple deprecated an algorithm that worked fine. From Maple 8, the algorithm moved to Rosenbrock methods for stiff equations. This is still not ideal.

If Maple had a working algorithm, I am hoping that Maplesoft folks would consider bringing it back in future versions. (At least with the same functionality as in Maple 6).

PLEASE NOTE, the issue is not with solving this example (Very simple). This example is chosen to show how a popular algorithm in the literature is wrongly implemented.

Here Maple's lsode is forced to take only one step and use first order back ward difference formula to integrate from 0 to 1. LSODE mimics Eulerbackward using the options given below. The post shows that LSODE does not do Newton Raphson and just performs only iteration for nonlinear equations.

>

restart;

>	Digits:=15;

(1)

>	eq:=diff(y(t),t)=-y(t);

(2)

>	C:=array([0$22]);

(3)

>

C[9]:=1;

(4)

>	sol:=dsolve({eq,y(0)=1},type=numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1):

>

sol(0.1);

(5)

>	subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq);

(6)

>	fsolve(%,y1=0.5);

(7)

While for linear it gave the expected result, it gives wrong results for nonlinear problems.

>	sol1:=dsolve({eq,y(0)=1},type=numeric):

>	sol1(0.1);

(8)

>	eq:=diff(y(t),t)=-y(t)^2exp(-y(t))-10y(t)(1+0.01exp(y(t)));

(9)

>	sol:=dsolve({eq,y(0)=1},type=numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1):

>

sol(0.1);

(10)

>	subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq);

(11)

>	fsolve(%,y1=1);

(12)

>	sol1:=dsolve({eq,y(0)=1},type=numeric):

the expected answer is correctly obtained with default tolerance as

>	sol1(0.1);

(13)

The results obtained are worse than single iteration using jacobian.

>	eq2:=(lhs-rhs)(subs(diff(y(t),t)=(y1-1)/0.1,y(t)=y1,eq));

(14)

>	jac:=unapply(diff(eq2,y1),y1);

(15)

>	f:=unapply(eq2,y1);

(16)

>

y0:=1;

(17)

>	dy:=-evalf(f(y0)/jac(y0));

(18)

>	ynew:=y0+dy;

(19)

Following procedures confirm that it is indeed calling the procedure only at 0 and 0.1, with backdiag giving slightly better results.

>	myfun:= proc(x,y) if not type(x,'numeric') or not type(evalf(y),numeric)then 'procname'(x,y); else lprint(`Request at x=`,x); -y^2exp(-y(x))-10y(1+0.01exp(y)); end if; end proc;

myfun := proc (x, y) if not (type(x, 'numeric') and type(evalf(y), numeric)) then ('procname')(x, y) else lprint(`Request at x=`, x); -y^2*exp(-y(x))-10*y*(1+0.1e-1*exp(y)) end if end proc

(20)

>	sol1:=dsolve({diff(y(x),x)=myfun(x,y(x)),y(0)=1},numeric,method=lsode[backfull],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1,known={myfun}):

>	sol1(0.1);

`Request at x=`, 0.

`Request at x=`, .1

(21)

>	sol2:=dsolve({diff(y(x),x)=myfun(x,y(x)),y(0)=1},numeric,method=lsode[backdiag],ctrl=C,initstep=0.1,minstep=0.1,abserr=1,relerr=1,known={myfun}):

>	sol2(0.1);

`Request at x=`, 0.

`Request at x=`, .1

(22)

Download Lsodeanalysistrunc.mws

Next see how dsolve method = mgear works just fine in Maple 6 (gives the expected answer upto 3 Digits accuracy). To run this code you will need Maple 6 or earlier versions. Maple 7 has this algorithm, but I don't know to use it as it is hidden. I would like to get support from other members to get Maplesoft's attention to bring this algorithm back.

If Mdy/dt = f(y) is solved using mgear algorithm (instead of dy/dt =f ), then one can have a good DAE solver based on this (M being singular).

Download Mgearworks.mws

Delay differential equations - Maple is ...

Posted: Dr. Venkat Subramanian 431

September 15 2015

8 6

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)^2exp(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));

eq3 := diff(y1(t), t) = piecewise(t <= .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)^2exp(1-y2(t)));

eq4 := diff(y2(t), t) = piecewise(t <= .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:

>

>	radar5data:=readdata("C:\\Users\\Venkat16core-office\\Google Drive\\waltmanproblem\\sol.txt",[string,string,float,string,string,float,float,float,float,float,float]):

>	nops(radar5data);

(7)

>	radar5data[1059];

(8)

>	eq[1]:=diff(y[1](t),t)=-ry[1](t)y[2](t)-sy[1](t)y[4](t);

(9)

>	eq[2]:=diff(y[2](t),t)=-ry[1](t)y[2](t)+alphary[1](y[5](t))y[2](y[5](t))H1;#Heaviside(t-35);

(10)

>	eq[3]:=diff(y[3](t),t)=ry[1](t)y[2](t);

(11)

>	eq[4]:=diff(y[4](t),t)=-sy[1](t)y[4](t)-gamma1y[4](t)+betary[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);

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

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

eq[6] := diff(y[6](t), t) = H2*(y[2](t)+y[3](t))/(y[2](y[6](t))+y[3](y[6](t)))

(14)

>	H1:=1/2+1/2tanh(100(t-35));H2:=1/2+1/2tanh(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;

[t = 300., y[1](t) = 0.615611371327094e-15, y[2](t) = 0.337706811581908e-6, y[3](t) = 0.422136411682798e-6, y[4](t) = 0.214253771204037e-5, y[5](t) = 299.999986716780, y[6](t) = 299.643054284209]

(19)

>	with(plots):

>	nd:=nops(radar5data);

(20)

>	radar5data[nd];

(21)

Values at t = 300 match with expected results.

>	pr[1]:=plot([seq([radar5data[i][3],log(radar5data[i][6])/log(10)],i=1..nd)],style=point,color=green):

>	p[1]:=odeplot(sol,[t,log(y[1](t))/log(10)],0..300,axes=boxed,thickness=3):

>	display({pr[1],p[1]});

>	pr[2]:=plot([seq([radar5data[i][3],log(radar5data[i][7])/log(10)],i=1..nd)],style=point,color=green):

>	p[2]:=odeplot(sol,[t,log(y[2](t))/log(10)],0..300,axes=boxed,thickness=3,numpoints=1000):

>	display({pr[2],p[2]});

>	pr[3]:=plot([seq([radar5data[i][3],log(radar5data[i][8])/log(10)],i=2..nd)],style=point,color=green):

>

>	p[3]:=odeplot(sol,[t,log(y[3](t))/log(10)],0..300,axes=boxed,thickness=3):

>	display({pr[3],p[3]});

>	pr[4]:=plot([seq([radar5data[i][3],log(radar5data[i][9])/log(10)],i=496..nd)],style=point,color=green,view=[197..300,-9..-5]):

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

>	pr[5]:=plot([seq([radar5data[i][3],radar5data[i][10]],i=1..nd)],style=point,color=green):

>	p[5]:=odeplot(sol,[t,y[5](t)],0..300,axes=boxed,thickness=3):

>	display({pr[5],p[5]});

>	pr[6]:=plot([seq([radar5data[i][3],radar5data[i][11]],i=1..nd)],style=point,color=green):

>	p[6]:=odeplot(sol,[t,y[6](t)],0..300,axes=boxed,thickness=3):

>	display({pr[6],p[6]});

Download delayimmunetopost.mws

Tricking Maple to solve Differential Alg...

Posted: Dr. Venkat Subramanian 431

July 08 2015

10 37

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

(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/2tanh(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;

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

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

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

(15)

>	eq2:=j1+j2=iapp;

(16)

>	j1:=io1(2(1-y1(t))exp((0.5F/R/T)(z1(t)-phi1))-2y1(t)exp((-0.5F/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};

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

(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/2tanh(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;

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;

(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/2tanh(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;

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

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

$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)-2y\|\|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)-2z\|\|i(t)+z\|\|(i-1)(t))/h^2-(1-y\|\|i(t)^2)(exp(-z\|\|i(t)));od:

>	eq1[1]:=(3y1(t)-4y2(t)+y3(t))/(2*h)=0: eq1[N+2]:=y\|\|(N+2)(t)-1=0:

>	eq2[1]:=(3z1(t)-4z2(t)+1z3(t))/(2h)=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/2tanh(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/2tanh(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

Download SolvingDAEsinMaple.mws

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

Venkat Subramanian

431 Reputation

13 Badges

MaplePrimes Activity

Maple helps in simulating multiscale-mul...

A new and efficient code in Maple for Bo...

Please bring back dsolve numeric method ...

Delay differential equations - Maple is ...

Tricking Maple to solve Differential Alg...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again