Items tagged with stiff

Feed

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;

Digits := 15

(1)

 

Enter the first-order ODEs into EqODEs list.

EqODEs:=[diff(y1(x),x)=y2(x),diff(y2(x),x)=epsilon*sinh(epsilon*y1(x))];

EqODEs := [diff(y1(x), x) = y2(x), diff(y2(x), x) = epsilon*sinh(epsilon*y1(x))]

(2)

 

Define the left boundary condition (bc1), and the right boundary condition (bc2). One should collect all the terms in one side.

bc1:=evalf([y1(x)]);

bc1 := [y1(x)]

(3)

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

bc2 := [y1(x)-1.]

(4)

 

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

Range:=[0.,1.];

Range := [0., 1.]

(5)

 

List any known parameters in the list.

pars:=[epsilon=2];

pars := [epsilon = 2]

(6)

 

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

unknownpars:=[];

unknownpars := []

(7)

 

Define the initial derivative in nder (default is 5 for 10th order) and the number of the nodes in nele (default is 10 and distributed evenly across the range provided by the user). The code adapts to increase the order. For many problems, 10th order method with 10 elements are sufficient.

nder:=5;nele:=10;

nder := 5

nele := 10

(8)

 

Define the absolute and relative tolerance for the local error. The error calculation is done based on the norm of both the 9th and 10th order simulation results.

atol:=1e-6;rtol:=atol/100;

atol := 0.1e-5

rtol := 0.100000000000000e-7

(9)

 

Call HDMadapt procedure, input all the information entered above and save the solution in sol. HDMadapt procedure does not need the initial guess for the mesh.

sol:= HDMadapt(EqODEs,bc1,bc2,pars,unknownpars,nder,nele,Range,atol,rtol):

 

Present some details of the solution.

sol[4]; # final derivative

5

(10)

sol[5]; # Maximum local RMSE

0.604570329905172e-8

(11)

 

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

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

NN := 11

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

node := 2

odevars := [y1(x), y2(x)]

(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

 

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!

 

``

deq1 := 1/(b-f(b)) = (2*(3-(1-f(b)*(diff(f(b), b, b)))/((diff(f(b), b))*(diff(f(b), b)))))/(1-2*(b-(1-f(b))/(diff(f(b), b))))

1/(b-f(b)) = 2*(3-(1-f(b)*(diff(diff(f(b), b), b)))/(diff(f(b), b))^2)/(1-2*b+2*(1-f(b))/(diff(f(b), b)))

(1)

ic1 := eval(f(b), b = 3/8) = 0, eval(f(b), b = 1/2) = 1/2

f(3/8) = 0, f(1/2) = 1/2

(2)

digits := 3

3

(3)

dsol1 := dsolve({deq1, ic1}, method = bvp[middefer], numeric, range = 3/8 .. 1/2)

Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

 

``

 

Download diffeqn.mw

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

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

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

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

epsilon := 0.1e-3

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

[t = 5.0, y1(t) = 1.60942323180838, y2(t) = .199998786891688]

1.609437912

.2000000000

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

1059

(7)

radar5data[1059];

["X", "=", 300.000000, "Y", "=", 0.6154486288e-15, 0.3377120916e-6, 0.4221403310e-6, 0.2142554563e-5, 299.9999999, 299.6430338]

(8)

eq[1]:=diff(y[1](t),t)=-r*y[1](t)*y[2](t)-s*y[1](t)*y[4](t);

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

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

(10)

eq[3]:=diff(y[3](t),t)=r*y[1](t)*y[2](t);

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

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

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

H1 := 1/2+(1/2)*tanh(100*t-3500)

H2 := 1/2+(1/2)*tanh(100*t-19700)

(15)

alpha:=1.8;beta:=20.;gamma1:=0.002;r:=5.*10^4;s:=10.^5;

alpha := 1.8

beta := 20.

gamma1 := 0.2e-2

r := 50000.

s := 100000.

(16)

seq(eq[i],i=1..6);

diff(y[1](t), t) = -50000.*y[1](t)*y[2](t)-100000.*y[1](t)*y[4](t), diff(y[2](t), t) = -50000.*y[1](t)*y[2](t)+90000.0*y[1](y[5](t))*y[2](y[5](t))*(1/2+(1/2)*tanh(100*t-3500)), diff(y[3](t), t) = 50000.*y[1](t)*y[2](t), diff(y[4](t), t) = -100000.*y[1](t)*y[4](t)-0.2e-2*y[4](t)+1000000.*y[1](y[6](t))*y[2](y[6](t))*(1/2+(1/2)*tanh(100*t-19700)), diff(y[5](t), t) = (1/2+(1/2)*tanh(100*t-3500))*(y[1](t)*y[2](t)+y[3](t))/(y[1](y[5](t))*y[2](y[5](t))+y[3](y[5](t))), diff(y[6](t), t) = (1/2+(1/2)*tanh(100*t-19700))*(y[2](t)+y[3](t))/(y[2](y[6](t))+y[3](y[6](t)))

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

ics := y[1](0) = 0.5000000000e-5, y[2](0) = 0.1000000000e-14, y[3](0) = 0, y[4](0) = 0, y[5](0) = 0.1e-39, y[6](0) = 0.1e-19

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

.141

(19)

with(plots):

nd:=nops(radar5data);

nd := 1059

(20)

radar5data[nd];

["X", "=", 300.000000, "Y", "=", 0.6154486288e-15, 0.3377120916e-6, 0.4221403310e-6, 0.2142554563e-5, 299.9999999, 299.6430338]

(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

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

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

xini := .1

.1

(1)

xfin := 2

2

(2)

SystemToSolve := diff(f[0, 0](x), x)+(2./x^5+.5000000000/x)*f[0, 0](x) = -15.58845727*sin(.5773502693*x)/x^2+(46.76537182*(3.*sin(.5773502693*x)-1.732050808*x*cos(.5773502693*x)))/x^4, diff(f[1, 0](x), x)+(6./x^5+1.500000000/x)*f[1, 0](x)-1.*f[0, 0](x)/x = (-15.58845727*sin(.5773502693*x)/x^2+(46.76537182*(3.*sin(.5773502693*x)-1.732050808*x*cos(.5773502693*x)))/x^4)*(1.-1.*(1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4))-(1.*(-10.39230485*sin(.5773502693*x)/x^2+(31.17691454*(3.*sin(.5773502693*x)-1.732050808*x*cos(.5773502693*x)))/x^4+(4.*((.8660254040*(3.*sin(.5773502693*x)-1.732050808*x*cos(.5773502693*x)))/x+(.8660254040*((3.*(1.-6./x^2))*sin(.5773502693*x)+10.39230485*cos(.5773502693*x)/x))/x))/((1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4)*x^5)))*(1/x^4)^(1/4)*exp(1/x^4)*GAMMA(.7500000000, 1/x^4), f[0, 0](xini) = 1.503498546, f[1, 0](xini) = -.5011661819:

 

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

 

``

ListProcs := dsolve({SystemToSolve}, numeric, method = lsode[backfull], output = listprocedure):

f00 := eval(f[0, 0](x), ListProcs);

proc (x) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x) else _xout := evalf(x) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := array( 1 .. 34, [( 1 ) = (2), ( 2 ) = (.1), ( 3 ) = (.1), ( 4 ) = (1), ( 5 ) = (1), ( 6 ) = (22), ( 7 ) = (0), ( 9 ) = (-.5011661819), ( 8 ) = (1.503498546), ( 11 ) = (0.1e-6), ( 10 ) = (0.1e-6), ( 13 ) = (0), ( 12 ) = (0), ( 15 ) = (0), ( 14 ) = (0), ( 18 ) = (0), ( 19 ) = (0), ( 16 ) = (0), ( 17 ) = (0), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = (0), ( 21 ) = (0), ( 27 ) = (0), ( 26 ) = (0), ( 25 ) = (0), ( 24 ) = (0), ( 31 ) = (0), ( 30 ) = (0), ( 29 ) = (0), ( 28 ) = (0), ( 32 ) = (0), ( 33 ) = (-1), ( 34 ) = (0)  ] ); _octl := array( 1 .. 34, [( 1 ) = (2), ( 2 ) = (.1), ( 3 ) = (.1), ( 4 ) = (1), ( 5 ) = (1), ( 6 ) = (22), ( 7 ) = (0), ( 9 ) = (-.5011661819), ( 8 ) = (1.503498546), ( 11 ) = (0.1e-6), ( 10 ) = (0.1e-6), ( 13 ) = (0), ( 12 ) = (0), ( 15 ) = (0), ( 14 ) = (0), ( 18 ) = (0), ( 19 ) = (0), ( 16 ) = (0), ( 17 ) = (0), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = (0), ( 21 ) = (0), ( 27 ) = (0), ( 26 ) = (0), ( 25 ) = (0), ( 24 ) = (0), ( 31 ) = (0), ( 30 ) = (0), ( 29 ) = (0), ( 28 ) = (0), ( 32 ) = (0), ( 33 ) = (-1), ( 34 ) = (0)  ] ); _n := trunc(_ctl[1]); _yini := Array(0..2, {(1) = .1, (2) = 1.503498546}); _y0 := Array(0..2, {(1) = .1, (2) = 1.503498546}); _fcn := proc (N, X, Y, YP) option `[Y[1] = f[0,0](x), Y[2] = f[1,0](x)]`; YP[1] := -15.58845727*sin(.5773502693*X)/X^2+46.76537182*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4-(2./X^5+.5000000000/X)*Y[1]; if 1/X^4 < 0 then YP[1] := undefined; return 0 end if; if 1/X^4 < 0 then YP[1] := undefined; return 0 end if; YP[2] := (-15.58845727*sin(.5773502693*X)/X^2+46.76537182*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4)*(1.-1.*evalf((1/X^4)^(1/4))*exp(1/X^4)*GAMMA(.7500000000, 1/X^4))-1.*(-10.39230485*sin(.5773502693*X)/X^2+31.17691454*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4+4.*(.8660254040*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X+.8660254040*(3.*(1.-6./X^2)*sin(.5773502693*X)+10.39230485*cos(.5773502693*X)/X)/X)*evalf(1/(1/X^4)^(1/4))/(exp(1/X^4)*GAMMA(.7500000000, 1/X^4)*X^5))*evalf((1/X^4)^(1/4))*exp(1/X^4)*GAMMA(.7500000000, 1/X^4)-(6./X^5+1.500000000/X)*Y[2]+1.*Y[1]/X; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "lsode" elif _xout = "numfun" then return trunc(_ctl[24+trunc(_ctl[1])]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[2]-_y0[0] = 0. then error "no information is available on last computed point" else _xout := _ctl[2] end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to _n do _octl[_i+7] := _y0[_i] end do; for _i to nops(_pars) do _octl[2*_n+30+_i] := _y0[_n+_i] end do; for _i to 34 do _ctl[_i] := _octl[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; for _i to _n do _octl[_i+7] := _y0[_i] end do; for _i to nops(_pars) do _octl[2*_n+30+_i] := _y0[_n+_i] end do end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[2] = 0 then [_ctl[2], seq(_ctl[_i], _i = 8 .. 7+_n)] else if sign(_xout-_ctl[2]) <> sign(_ctl[2]-_y0[0]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 34 do _ctl[_i] := _octl[_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/lsode`(_fcn, var(_ctl))) catch: userinfo(2, `dsolve/debug`, print(`Exception in lsode:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/lsode`(_fcn, _ctl) else error  end if end try else _errcd := `dsolve/numeric/lsode`(_fcn, _ctl) end if; if _errcd < 0 then userinfo(2, {dsolve, `dsolve/lsode`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/lsode`}, ` t =`, _ctl[2]); _i := 8; userinfo(2, {dsolve, `dsolve/lsode`}, ` y =`, _ctl[_i]); for _i from _i+1 to 7+_n do userinfo(2, {dsolve, `dsolve/lsode`}, `	 `, _ctl[_i]) end do; if _errcd+1. = 0. then if _ctl[14+trunc(_ctl[1])] <> 0 then error "an excessive amount of work was done, maxstep may be too small" else error "an excessive amount of work (greater than mxstep) was done" end if elif _errcd+2. = 0. then error "too much accuracy was requested for the machine being used" elif _errcd+3. = 0. then error "illegal input was detected" elif _errcd+4. = 0. then error "repeated error test failures on the attempted step" elif _errcd+5. = 0. then error "repeated convergence test failures on the attempted step" elif _errcd+6. = 0. then error "pure relative error control requested for a variable that has vanished" elif _errcd+7. = 0. then error "cannot evaluate the solution past %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_ctl[2]) else error "unknown error code returned from lsode %1", trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_ctl[_i], _i = 8 .. 7+_n)] end if end proc, (2) = Array(1..3, {(1) = 18446744078356217278, (2) = 18446744078356217454, (3) = 18446744078356217630}), (3) = [x, f[0, 0](x), f[1, 0](x)], (4) = []}); _solnproc := _dat[1]; _pars := map(rhs, _dat[4]); if not type(_xout, 'numeric') then if member(x, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x, ["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(x, 'string')); if type(_res, 'list') then return _res[2] else return NULL end if elif member(x, ["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x), 'string') = rhs(x); if lhs(_xout) = "initial" then if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 2, rhs(_xout)]) end if elif not type(rhs(_xout), 'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 2, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[2] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[2], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x), 'string') = rhs(x)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(x) else _ndsol := `tools/gensym`("f[0,0](x)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"), _Inert_EXPSEQ(ToInert(_ndsol), _Inert_VERBATIM(pointto(_dat[2][2])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol), _Inert_EXPSEQ(ToInert(x)))) end if end if; try _res := _solnproc(_xout); _res[2] catch: error  end try end proc

(3)

f10 := eval(f[1, 0](x), ListProcs);

proc (x) local _res, _dat, _solnproc, _xout, _ndsol, _pars, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x) else _xout := evalf(x) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := array( 1 .. 34, [( 1 ) = (2), ( 2 ) = (.1), ( 3 ) = (.1), ( 4 ) = (1), ( 5 ) = (1), ( 6 ) = (22), ( 7 ) = (0), ( 9 ) = (-.5011661819), ( 8 ) = (1.503498546), ( 11 ) = (0.1e-6), ( 10 ) = (0.1e-6), ( 13 ) = (0), ( 12 ) = (0), ( 15 ) = (0), ( 14 ) = (0), ( 18 ) = (0), ( 19 ) = (0), ( 16 ) = (0), ( 17 ) = (0), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = (0), ( 21 ) = (0), ( 27 ) = (0), ( 26 ) = (0), ( 25 ) = (0), ( 24 ) = (0), ( 31 ) = (0), ( 30 ) = (0), ( 29 ) = (0), ( 28 ) = (0), ( 32 ) = (0), ( 33 ) = (-1), ( 34 ) = (0)  ] ); _octl := array( 1 .. 34, [( 1 ) = (2), ( 2 ) = (.1), ( 3 ) = (.1), ( 4 ) = (1), ( 5 ) = (1), ( 6 ) = (22), ( 7 ) = (0), ( 9 ) = (-.5011661819), ( 8 ) = (1.503498546), ( 11 ) = (0.1e-6), ( 10 ) = (0.1e-6), ( 13 ) = (0), ( 12 ) = (0), ( 15 ) = (0), ( 14 ) = (0), ( 18 ) = (0), ( 19 ) = (0), ( 16 ) = (0), ( 17 ) = (0), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = (0), ( 21 ) = (0), ( 27 ) = (0), ( 26 ) = (0), ( 25 ) = (0), ( 24 ) = (0), ( 31 ) = (0), ( 30 ) = (0), ( 29 ) = (0), ( 28 ) = (0), ( 32 ) = (0), ( 33 ) = (-1), ( 34 ) = (0)  ] ); _n := trunc(_ctl[1]); _yini := Array(0..2, {(1) = .1, (2) = 1.503498546}); _y0 := Array(0..2, {(1) = .1, (2) = 1.503498546}); _fcn := proc (N, X, Y, YP) option `[Y[1] = f[0,0](x), Y[2] = f[1,0](x)]`; YP[1] := -15.58845727*sin(.5773502693*X)/X^2+46.76537182*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4-(2./X^5+.5000000000/X)*Y[1]; if 1/X^4 < 0 then YP[1] := undefined; return 0 end if; if 1/X^4 < 0 then YP[1] := undefined; return 0 end if; YP[2] := (-15.58845727*sin(.5773502693*X)/X^2+46.76537182*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4)*(1.-1.*evalf((1/X^4)^(1/4))*exp(1/X^4)*GAMMA(.7500000000, 1/X^4))-1.*(-10.39230485*sin(.5773502693*X)/X^2+31.17691454*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X^4+4.*(.8660254040*(3.*sin(.5773502693*X)-1.732050808*X*cos(.5773502693*X))/X+.8660254040*(3.*(1.-6./X^2)*sin(.5773502693*X)+10.39230485*cos(.5773502693*X)/X)/X)*evalf(1/(1/X^4)^(1/4))/(exp(1/X^4)*GAMMA(.7500000000, 1/X^4)*X^5))*evalf((1/X^4)^(1/4))*exp(1/X^4)*GAMMA(.7500000000, 1/X^4)-(6./X^5+1.500000000/X)*Y[2]+1.*Y[1]/X; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "lsode" elif _xout = "numfun" then return trunc(_ctl[24+trunc(_ctl[1])]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[2]-_y0[0] = 0. then error "no information is available on last computed point" else _xout := _ctl[2] end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to _n do _octl[_i+7] := _y0[_i] end do; for _i to nops(_pars) do _octl[2*_n+30+_i] := _y0[_n+_i] end do; for _i to 34 do _ctl[_i] := _octl[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; for _i to _n do _octl[_i+7] := _y0[_i] end do; for _i to nops(_pars) do _octl[2*_n+30+_i] := _y0[_n+_i] end do end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[2] = 0 then [_ctl[2], seq(_ctl[_i], _i = 8 .. 7+_n)] else if sign(_xout-_ctl[2]) <> sign(_ctl[2]-_y0[0]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 34 do _ctl[_i] := _octl[_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/lsode`(_fcn, var(_ctl))) catch: userinfo(2, `dsolve/debug`, print(`Exception in lsode:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/lsode`(_fcn, _ctl) else error  end if end try else _errcd := `dsolve/numeric/lsode`(_fcn, _ctl) end if; if _errcd < 0 then userinfo(2, {dsolve, `dsolve/lsode`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/lsode`}, ` t =`, _ctl[2]); _i := 8; userinfo(2, {dsolve, `dsolve/lsode`}, ` y =`, _ctl[_i]); for _i from _i+1 to 7+_n do userinfo(2, {dsolve, `dsolve/lsode`}, `	 `, _ctl[_i]) end do; if _errcd+1. = 0. then if _ctl[14+trunc(_ctl[1])] <> 0 then error "an excessive amount of work was done, maxstep may be too small" else error "an excessive amount of work (greater than mxstep) was done" end if elif _errcd+2. = 0. then error "too much accuracy was requested for the machine being used" elif _errcd+3. = 0. then error "illegal input was detected" elif _errcd+4. = 0. then error "repeated error test failures on the attempted step" elif _errcd+5. = 0. then error "repeated convergence test failures on the attempted step" elif _errcd+6. = 0. then error "pure relative error control requested for a variable that has vanished" elif _errcd+7. = 0. then error "cannot evaluate the solution past %1, maxfun limit exceeded (see <a href='http://www.maplesoft.com/support/help/search.aspx?term=dsolve,maxfun' target='_new'>?dsolve,maxfun</a> for details)", evalf[8](_ctl[2]) else error "unknown error code returned from lsode %1", trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_ctl[_i], _i = 8 .. 7+_n)] end if end proc, (2) = Array(1..3, {(1) = 18446744078356217278, (2) = 18446744078356217454, (3) = 18446744078356217630}), (3) = [x, f[0, 0](x), f[1, 0](x)], (4) = []}); _solnproc := _dat[1]; _pars := map(rhs, _dat[4]); if not type(_xout, 'numeric') then if member(x, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x, ["last", 'last', "initial", 'initial', NULL]) then _res := _solnproc(convert(x, 'string')); if type(_res, 'list') then return _res[3] else return NULL end if elif member(x, ["parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[3], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x), 'string') = rhs(x); if lhs(_xout) = "initial" then if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else _res := _solnproc("initial" = ["single", 3, rhs(_xout)]) end if elif not type(rhs(_xout), 'list') then error "initial and/or parameter values must be specified in a list" elif lhs(_xout) = "initial_and_parameters" and nops(rhs(_xout)) = nops(_pars)+1 then _res := _solnproc(lhs(_xout) = ["single", 3, op(rhs(_xout))]) else _res := _solnproc(_xout) end if; if lhs(_xout) = "initial" then return _res[3] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [_res[3], seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x), 'string') = rhs(x)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _dat[3] end if; if procname <> unknown then return ('procname')(x) else _ndsol := `tools/gensym`("f[1,0](x)"); eval(FromInert(_Inert_FUNCTION(_Inert_NAME("assign"), _Inert_EXPSEQ(ToInert(_ndsol), _Inert_VERBATIM(pointto(_dat[2][3])))))); return FromInert(_Inert_FUNCTION(ToInert(_ndsol), _Inert_EXPSEQ(ToInert(x)))) end if end if; try _res := _solnproc(_xout); _res[3] catch: error  end try end proc

(4)

ftoplot := unapply(f00(x)+0.45e-1*f10(x)/x^(3/2), x)

proc (x) options operator, arrow; f00(x)+0.45e-1*f10(x)/x^(3/2) end proc

(5)

``

plot(ftoplot(x), x = xini .. xfin)

 

``


Download Test2ODEs.mw

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.  

 

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:

 

Sollsodefull:=dsolve({ODEs}, numeric, method = lsode[adamsfull])

 

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

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?

I really appreciate your help.

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}:

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.

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

Dear mapleprimes ODE enthusiasts,

I am analyzing a stiff 3D system characterizes by fast-slow dynamics. I believe I can show analytically that, for some parameter values, in a neighborhood of some critical points (where the system is stationary), the dynamic system has a local center manifold that is center-stable. I have simulated the system for different parameter values.I am interested in one particular trajectory. With a random set of initial values, the system...

Page 1 of 1