Items tagged with dsolve dsolve Tagged Items Feed

dsolve thousands of system with for loop

but one of system run a very long time, seems never end, just stuck here

how to set timeout when impossible to dsolve

ics:=a(1)=0,b(1)=0,c(1)=0:
ode:=rhs(a4[ii])[1]=1,rhs(a4[ii])[2]=1,rhs(a4[ii])[3]=1:
ode:=a(t)*c(t)*(diff(c(t), t))+2*a(t)*(diff(a(t), t))*b(t) = 1, a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))*c(t)*(diff(c(t), t))+a(t)*(diff(a(t), t))*b(t)*(diff(b(t), t))+a(t)*(diff(a(t), t)) = 1, a(t)*(diff(a(t), t))*c(t)+2*a(t)*b(t)*(diff(b(t), t)) = 1;
dsol := dsolve({ode,ics}):

 

 

found differentiated functions with same name but depending on different arguments in the given DE system

restart:

eqn0 := diff(f0(y), y, y, y, y) = 0;

bcs0:=f0(0)=0,(D@@2)(f0)(0)=0,f0(S)=1,D(f0)(S)=0:

dsolve({eqn0,bcs0});

How to fix it?

Thanks

Hi,
I'm new at maple and have a problem/question with the rkf45 numerical ODE Solver.

At first, my computer need a lot of time to calculate an analytic solution.
Therefor, I use the numerical way.

I have the following second order ODE:
ODE:=m*((D@@2)(x))(t)+d*(D(x))(t)+k*x(t) = d*(eval(diff(y(x), x), x = t))+k*y(t)
where y(t) is a realy big piecewise function, defined by me.

My initial conditions are:
x(0) = 0, (D(x))(0) = 0

With dsolve, I get the solution x(t) and the first derivative x'(t). I'm able to plot them with odeplot.

But...

Problem 1:
I need also the second derivative x''(t).
On this page: http://www.maplesoft.com/support/help/maple/view.aspx?path=dsolve%2Frkf45there is an example (eq 13 and 14) where the second derivative is useable, but this doesn't work with my differential equation.
I have add 
(D(D(x)))(0) = 0
to my initial conditions but then, I got the error that only 2 initial conditions are required.
What could I do, so that rkf45 returns also the second derivative?

Problem 2:
And in addition to this, I want to calculate with x(t), x'(t), x''(t) but I found no way to use them.
Only plots are possible.
If I reduce y(t) to a minimum, I can do everything with the analytic solution: plot, d/dt, d2/dt2, +, -, ...
I tried also to convert the procedure to a function but in this case, there is no way to derivate it.

Many thanks...

I recently got myself a 64 bit computer and have noticed that I cannot use the option compile=true in dsolve/numeric. Take the following simple example:

dsolve({diff(x(t),t)=x(t),x(0)=1},numeric,compile=true);
Error, (in dsolve/numeric/SC/preproc) unable to compile (rc=1), please try again, and if that fails verify your Windows compiler installation

I'm using Windows 10, but had the same problem with Windows 8.1 on the same machine.
The Compiler:-Compile examples in the help page all work.

What do I have to do to make the option compile=true work in dsolve/numeric?

You may safely assume that I don't know any technicalities about these things.

restart;
odes:=diff(x1(t),t)*diff(x2(t),t$2)*sin(x1(t)*x2(t))+5*diff(x1(t),t$2)*diff(x2(t),t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2=exp(-x2(t)^2),diff(x1(t),t$2)*x2(t)+diff(x2(t),t$2)*diff(x1(t),t)*sin(x1(t)^2)+cos(diff(x2(t),t$2)*x2(t))=sin(t);
ics:=x1(0)=1,D(x1)(0)=1,x2(0)=2,D(x2)(0)=2;
subs(diff(x1(t),t$2)=yp2,diff(x2(t),t$2)=yp4,diff(x1(t),t)=Y[2],diff(x2(t),t)=Y[4],x1(t)=Y[1],x2(t)=Y[3],{odes});
p:=proc(N,t,Y,YP)
local eqs,yp2,yp4;
YP[1]=Y[2];
YP[3]=Y[4];
eqs:=[yp2*Y[3]+yp4*Y[2]*sin(Y[1]^2)+cos(yp4*Y[3]) = sin(t), Y[2]*yp4*sin(Y[1]*Y[3])+5*yp2*Y[4]*cos(Y[1]^2)+t^2*Y[1]*Y[3]^2 = exp(-Y[3]^2)];
YP[2],YP[4]:=op(subs(fsolve(eqs,{yp2=1,yp4=2}),[yp2,yp4]));
end proc:
res:=dsolve(numeric,procedure=p,initial=array([1,1,2,2]),number=4,procvars=[x1(t),diff(x1(t),t),x2(t),diff(x2(t),t)],start=0,maxfun=0);
plots:-odeplot(res,[t,x1(t)],0..5,gridlines=true);

It's a long time to wait for the odeplot.

Any advice is appreciated.

Hi all,

restart;#part1
epsilon:=5:Delta1:=2:Delta2:=-4: N1:=1000:
dsys :={diff(x(t),t)=-I*Delta1*x(t)+y(t)+epsilon, diff(y(t),t)=-I*Delta2*y(t)+x(t)*z(t), diff(z(t),t)=-2*(conjugate(x(t))*y(t)+conjugate(y(t))*x(t))};

res:=dsolve(dsys union {x(0)=2*I,y(0)=0,z(0)=1},numeric,output=listprocedure);

P1:=plots:-odeplot(res,[[t,Im(y(t))],[t,Re(x(t))]],0..2):

/ d
{ --- x(t) = -2 I x(t) + y(t) + 5,
\ dt

d
--- y(t) = 4 I y(t) + x(t) z(t),
dt

d ____ ____ \
--- z(t) = -2 x(t) y(t) - 2 y(t) x(t) }
dt /
tit:=sprintf("D1=%g,D2=%g",Delta1,Delta2);
"D1=2,D2=-4"
plots[odeplot](res,[[t,Im(y(t))]],0..200,axes=boxed,titlefont=[SYMBOL,14],font=[1,1,18],color=black,linestyle=1,tickmarks=[3, 4],font=[1,1,14],thickness=2,titlefont=[SYMBOL,12]);
Warning, cannot evaluate the solution further right of 90.013890, maxfun limit exceeded (see ?dsolve,maxfun for details)

when I increase the time give this msn:

Warning, cannot evaluate the solution further right of 90.013890, maxfun limit exceeded (see ?dsolve,maxfun for details)

hi.I want to dsolve set of nonlinear equations with one unknown parameter ...is this possible with dsolve rule.in matlab this possible with bvp4c rule..please help me for this problem.if we should another rule please attached file reform.Thanks alot12.mw

Hi,

I am trying to solve a set of ode which depends on some parameters like A0,0,  A1,0, A1,1, B0,0,  B1,0, B1,1, C0,0 and so on. Here is some part of my code:

 

restart;

sigma := 1; X := proc (m) if m <= 1 then 0 elif 2 <= m then 1 end if end proc;

lambda := proc (m, k) if k <= m and 0 <= k then 1 else 0 end if end proc;

Typesetting:-Settings(functionassign = false);

for m from 0  to 4 do    

for k  from 4 by -1 to 0 do  

A[m,k](r) :=diff(f[m,k](r),r$2)+((k+1)*(k+2))/(r^(2))*(lambda(m,k+2)*f[m,k+2](r)-f[m,k](r)):  

T[m,k]:=(k+1)*(lambda(m,k+1)*A[m,k+1](r) -lambda(m,k-1)*A[m,k-1](r)):      

S[m,k]:=(k+1)*(lambda(m,k+1)*f[m,k+1]( r)-lambda(m,k-1)*f[m,k-1](r)):               # There are also C[m,k], B[m,k] and E[m,k] definitions similar to A[m,k],T[m,k] and S[m,k]

Q[m, k] := r^(4-sigma)*E[m, k]-2*lambda(m, k+2)*(k+1)*(k+2)*(r^2*(diff(f[m, k+2](r), `$`(r, 2)))-2*r*(diff(f[m, k+2](r), r)))-(k+1)*(k+2)*(k+4)*((k+3)*lambda(m, k+4)*f[m, k+4](r)-(2*(k+1))*lambda(m, k+2)*f[m, k+2](r));

ode[m, k] := r^4*(diff(f[m, k](r), `$`(r, 4)))-(k+1)*(k+2)*(2*r^2*(diff(f[m, k](r), `$`(r, 2)))-4*r*(diff(f[m, k](r), r))-(k-1)*(k+4)*f[m, k](r)) = Q[m, k];

soln[m, k] := rhs(dsolve(ode[m, k], f[m, k](r)))

 

After obtaining the solution, there are some additional parts in my code to find the coefficients of odes. The code is just working fine until m=3. When m=3, I am getting this error:

Error, (in solve) cannot solve expressions with Int(-(1/282240)*r*(3*(h . R)^2*(Int((6720*ln(r)*h*r^6-22176*ln(r)*h*r^5-50400*ln(r)*r^6-2814*h*r^6+33600*r^7+18144*ln(r)*h*r^4+90720*ln(r)*r^5+8238*h*r^5-95520*r^6+2352*h*ln(r)*r^3-7974*h*r^4+57780*r^5-2880*ln(r)*h*r^2-14400*ln(r)*r^3-2443*h*r^3+2100*r^4+6108*h*r^2+5340*r^3-1115*h-3300*r)/r^6, r))-4480*_C1), r) for _C1

I think the problem is due to declaration of f[m,k](r) values for the set of odes. For example, before solving ode[3,3], the code declares a value for f[3,3](r) which includes some integral definition in ode[3,3] although i want to find the solution for f[3,3](r). 

To illustrate,

ode[1,1]=r^4*(diff(f[1, 1](r), r, r, r, r))-12*r^2*(diff(f[1, 1](r), r, r))+24*r*(diff(f[1, 1](r), r)) = -Typesetting[delayDotProduct](r, h . R, true)*((2*(diff(f[0, 0](r), r))-4*f[0, 0](r)/r)*(diff(f[0, 0](r), r, r)-2*f[0, 0](r)/r^2)-(2*(-3/4-1/(4*r^2)+r))*A[0, 0]+(2*(diff(f[0, 0](r), r, r, r)+4*f[0, 0](r)/r^3-2*(diff(f[0, 0](r), r))/r^2))*f[0, 0](r))

soln[m, k] := rhs(dsolve(ode[m, k], f[m, k](r)))

 

where I already know the definition of f[0,0](r)=-(3/4)*r+1/(4*r)+(1/2)*r^2   and A[0,0]=1/(2*r^3)+1+(2*((3/4)*r-1/(4*r)-(1/2)*r^2))/r^2 

So, ode[1,1] can be solved with respect to f[1,1](r).

I would be glad for any comments. 

Here is my code.

 test.mw

after tried dsolve(sol, [a,b,c]) or dsolve(sol,t), still have errors

First equation

sol := a*(diff(M(a, b, c), a))+a*b*(diff(M(a, b, c), b))*c
dsolve(sol);
Error, (in ODEtools/info) Required a specification of the indeterminate function

Second Equation

sol := a(t)*(diff(a(t), t))+a(t)*b(t)*(diff(b(t), t))*c(t)
dsolve(sol);
Error, (in ODEtools/info) Required a specification of the indeterminate function

Hi !

I wish to solve a differential equation with an initial condition of the type D(f)(a)=b, where 'a' and 'b' are arbitrary real constants and are entered this way symbolically. Maple interprets D(f)(a) as a function which depends on a new variable 'a' and gives the error:

 

Error, (in dsolve) found differentiated functions with same name but depending on different arguments in the given DE system: {f(a), f(psi)}

(the differential equation is for a function f(psi) which was defined previously)

 

How can I pass this initial condition f'(a)=b (a,b being constants, but arbitrary) to the Maple dsolve command?

I have been browsing the web and looking at manuals, but always find D(f)(0)=3 or any other combination with numbers only.

Thank you in advance.

Regards,

              Michael

 

 

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 make an attempt to plot the solution to

Here is my code :

> with(plots); with(DEtools);
> ode1 := diff(x(t), t) = v(t); ode2 := diff(v(t), t) = -(.8*9.8)*v(t)/abs(v(t))-cos(t)^2;
> MODEL := {ode1, ode2}; VARS := {v(t), x(t)}; DOMAIN := t = 0 .. 150; RANGE := x = -1 .. 1, v = -5 .. 5; COLORS := [BLACK, BLUE]; IC1 := [x(0) = .5, v(0) = .25]; IC2 := [x(0) = 2.5, v(0) = 3];
> DEplot(MODEL, VARS, DOMAIN, RANGE, [IC1, IC2], stepsize = .1, linecolor = COLORS, scene = [t, x]);
>

and the message cannot evaluate the solution further right of .16015784, maxfun limit exceeded (see ?dsolve,maxfun for details)

Any other attemp has failed.

Have you got somme ideas

Thanks

Phil

Hello, dear experts.
I have a question...
solve the system of differential equations,where one of the initial conditions need to be chosen so thatcondition is metat the end of integration.
The task is not difficult, but I'm having trouble with the syntax.

1.I can't "pull"the desired function from the solution and find its value at a certain point.
I try to do so:
r_ravn:=s->subs(F,r(s));
evalf(r_ravn(s_end));
evalf(r_ravn(0));
but there is no result

2.In this case,instead of"for"it is better to use a while loop, but again the problem arises 1.
Tell me, please,how to implemen my program.

 

restart:
R:=0.3:
theta_min:=Pi/6:
theta_max:=Pi/2:
betta_max:=evalf(Pi/180*80);
p:=2*10^5:

theta0:=s->Pi/3/s_end*s+Pi/6:
r0:=s->R*sin(theta0(s)):
s_end:=evalf(R*(theta_max-theta_min)):

sol1:=solve({sin(betta_max)=c/r0(0)},{c});
const1:=0.1477211630;

betta0:=s->arcsin(const1/r0(s)):
betta:=s->arcsin(r(s)/r0(s)*sin(betta0(s))):
A:=s->cos(betta(s))/cos(betta0(s)):
T1:=s->rT1(s)/r(s):
T2:=s->T1(s)*tan(betta(s))^2:

step:=0.001:
delta:=0.001:
for i from 1 to 3000 do
r_min:=0.3-step:
rT1_n:=p*Pi*r_min^2/2/Pi/sin(theta_min):

sys := diff(rT1(s), s)-A(s)*T2(s)*cos(theta(s)),diff(theta(s), s)-A(s)/T1(s)*(p-T2(s)*sin(theta(s)/r(s))),diff(r(s),s)-A(s)*cos(theta(s)),diff(z(s),s)-A(s)*sin(theta(s));
fcns := {rT1(s),theta(s),r(s), z(s)};
F := dsolve({sys,rT1(0)=rT1_n, theta(0)=theta_min,r(0) = r_min, z(0) = 0}, fcns, numeric,output=listprocedure):
r_ravn:=s->subs(F,r(s)):
if abs(evalf(r_ravn(s_end))-R)=delta then break:
print(r_min):
end if:
end do:

r_ravn:=s->subs(F,r(s));
evalf(r_ravn(s_end));
evalf(r_ravn(0));
plot([r_ravn(s),r(s)],s=0..s_end);

Here my worksheet that at the end of it I have a problem with the dsolve when solving an ode system. and the error is 

Error, (in StringTools:-IsPrefix) second argument must be a string

I know that the dsolve has problem with bracket . Then how can I fix it or how can I change my codes.

 optimal.mwoptimal.mw

I'm currently having some difficulties in solving a system of differential equations numerically.

This is my code.

 

1 2 3 4 5 6 7 Last Page 1 of 23