Advanced Search

Venkat Subramanian

411 Reputation

13 Badges

14 years, 343 days
Contact 119368

MaplePrimes Activity


These are answers submitted by

good discussion and alternate approach...

Dr. Venkat Subramanian 411
August 25 2023
1 0

30 minutes ago

@Rouben Rostamian  , @rlopez

Great discussion.
For @kjell
Please see codes at 
https://sites.utexas.edu/maple/models-codes-2/challenges-in-moving-to-multiscale-battery-models-where-electrochemistry-meets-and-demands-more-from-math/ 
The codes involved biquadratic FEM for Laplace equation with singular corners with rectangular elements.
As you might know for large number of elements, you will need a direct sparse linear solver (robust) or iterative Krylov type solvers(fast, but needs good preconditioner and may not be robust). The codes used SparseLU. The newer Maple 2023 is even better and has SpasreMKLDirect for intel processors.

An alternative approach to solve irregular domains is the method of immersed domain. That is you wrap irregular domains and two phases into one rectangle and then solve.
In this case, Laplace equation becomes Poisson's equation (addition of a source term). Typically this approach is only approximate at the interfaces and is quite reasonable away from the interfaces. This forms the basis for moving boundary models, phase-field models, etc.
See the discussion at https://www.mapleprimes.com/questions/234050-Will-Anyone-Be-Able-To-Speed-Up-This 
With the help of @acer, I was able to optimize calls to UMFPACK and solve 200x200 grid + time stepping in few seconds.
Later, with the help of @epostma, I was able to use SparseMKLDirect to solve even 1000x1000 grid with 100 + time steps in few seconds or less including cycling of batteries. Paper and code at 
https://sites.utexas.edu/maple/phase-field-models-for-li-metal-batteries/ 

We also developed a Maple solver for largescale DAE problems (under review at Maple Transactions). Paper and codes at 
https://sites.utexas.edu/maple/sparsedaesolver/ 

This solver can easily solve upto 10000 DAEs. For 1M DAEs, you have to use the MKL sparse solvers (available with 2023). If you want you can add a dummy time derivative to your equation and take it to steadystate.


Time permitting I plan to provide interfaces to triangular meshing and FEM packages as well. Note that for general PDEs, it is not just triangulation, one needs to write the residues and discretized equations. Maple can play a big role in the same.

Maple is indeed very powerful. IMHO, it is the "easiest to code, and closest to C for efficiency".

Another (perhaps faster) approach -Simul...

Dr. Venkat Subramanian 411
April 11 2021
0 0

The original problem has only 3 equations. I rewrote the model to write 3 equations for every value of delta (I used 21 data points for delta, you can increase it, but it might takDownload SimultaneousSolution.mwe more RAM).

Hence solving for 3*21 = 63 equations (21 different values of delta) seems to give faster results. Following CPU time is reported (includes time taken to dsolve and one plot. You should be able to return data till t= 0.5 and make 3D plots).

memory used=0.82GiB, alloc change=0 bytes, cpu time=5.84s, real time=5.85s, gc time=1.09s
 

If you change the model to real part and imaginary part, I expect dsolve to be much faster and RAM efficient.

 

PS., if anyone edits this, please avoid using ~/, shortcuts so that code is easy to follow and reverse compatible with older versions of Maple. (I am willing to learn shortcuts if gives me gain in cpu time/efficiency/RAM usage).

I am not able to display contents from my code. So, copying and pasting relevant lines of codes below.

restart;
Digits:=15;
                               15

phi:=0:lambda:=0.1:N:=5:M:=sqrt(N*(N+1))*exp(I*phi):omegap:=10:
dsys:=[diff(n(t),t)=-2*(n(t)-N)+(u(t)-M)*exp(-2*I*omegap*t/lambda)
       +((z(t)-conjugate(M))*exp(2*I*omegap*t/lambda)),
       diff(u(t),t)=-2*(1-I*delta)*u(t)
                    +2*(n(t)-N)*exp(2*I*omegap*t/lambda)+2*M,
       diff((z(t),t))=-2*(1+I*delta)*z(t)
                      +2*(n(t)-N)*exp(-2*I*omegap*t/lambda)
                      +2*conjugate(M)];
 deltavals:=[seq(i-10.,i=0..20)];
[-10., -9., -8., -7., -6., -5., -4., -3., -2., -1., 0., 1., 2., 

  3., 4., 5., 6., 7., 8., 9., 10.]


Ntot:=nops(deltavals);
                               21

vars:=[n(t)=NN[i1](t),u(t)=UU[i1](t),z(t)=ZZ[i1](t)];
     [n(t) = NN[i1](t), u(t) = UU[i1](t), z(t) = ZZ[i1](t)]

for i from 1 to Ntot do eqs[i]:=evalf(subs(vars,i1=i,delta=deltavals[i],dsys));od:
ics:=seq(eval(NN[j](0))=0,j=1..Ntot),seq(eval(UU[j](0))=0,j=1..Ntot),seq(eval(ZZ[j](0))=0,j=1..Ntot):


Dsystot:=seq(op(eqs[i]),i=1..Ntot):

soln:=dsolve({Dsystot,ics},numeric,compile=true):
ss:=soln(0.5):

plot([seq([deltavals[i],subs(ss,NN[i](t))],i=1..Ntot)],style=point);


A simple procedure is written to find the cpu time/RAM usage for find the dsolve solution and then creating data for plots. Once sol(0.5) is found, solution at t= 0.5,
Simul:=proc()
local soln1,ss1,i;
soln1:=dsolve({Dsystot,ics},numeric,compile=true):
ss1:=soln1(0.5):
plot([seq([deltavals[i],subs(ss,NN[i](t))],i=1..Ntot)],style=point);
end proc;

with(CodeTools:-Profiling):
CodeTools:-Usage(Simul());
memory used=0.82GiB, alloc change=0 bytes, cpu time=5.84s, real time=5.85s, gc time=1.09s

 

 

see this post...

Dr. Venkat Subramanian 411
June 27 2020
2 5

and use a piecewise function

 

https://www.mapleprimes.com/posts/201287-Delay-Differential-Equations--Maple-Is-A-Leader

 

still slow...

Dr. Venkat Subramanian 411
May 15 2020
0 6

@janhardo 

If you are used to old mws format, then the best you can do today is to use worksheet mode with 1D input.

Note that current worksheet mode (from Maple 10 I believe) is Java based, and is very slow compared to classic worksheet model.

You can run current 32 bit versions in classic worksheet mode, but it crashes very frequently deleting all the commands. Dell seems to be better than hp computers for stability (in windows 10). Microsoft tablet is the worst (I could never install any stable mws format in that) and threw that laptop away.

fsolve sucks big time...

Dr. Venkat Subramanian 411
March 05 2020
1 3

If you have a good initial guess, try to write a simple Newton Raphson code. 

As of now, fsolve will try to increase Digits arbitrarily and won't solve for nonlinear equations to a  specified tolerance. fsolve works only for nearly perfect residuals (= 0). For a larger system, it may not give any answers.

95 equations should be solvable. You can also try Optimization:-NLPSolve(1, []) with equations as constraints.

 

Define a new and equivalent BVP such tha...

Dr. Venkat Subramanian 411
February 13 2018
1 1

I don't recommend Shooting methods. In particular, single shooting method will fail for many nonlinear and stiff problems. Multiple shooting is better. Use of optimization (combining IVP or BVP solver with an optimizer) is inefficient and should be used only for optimal control or parameter estimation unless you know that the BVPs have a reasonable initial guess for unknown boundary conditions at one end.

Divide x = 0.. 2 

as x = 0..1 (region 1) , and x2 = x-1 = 0..1 (region 2)

Define y as y in region 1 and y2 in region 2. 

Without loss of generality replace x2 with x. Make sure you don't have any explicit functions of x (non-homogeneous non-constant terms in your equation. If there is any, use a dummy variable z and use dz/dx = 1)

We need additional boundary conditions. Derivatives and second derivatives are continuous at region 1/region 2 interface. So

dy/dx(1) = dy2/dx(0), etc.

For details on related approach (to scale any domain to 0 to 1), see past papers on BVP solvers or one of our past papers, (figure 1 in particular)

http://depts.washington.edu/maple/pubs/40_JES_OC_reformulation.pdf

You might want to see a different approach to solve BVPs (in particular if Maple's BVP solver fails or if you want to solve a largescale BVP).

https://www.mapleprimes.com/posts/208499-A-New-And-Efficient-Code-In-Maple-For

 

> 

restart;
> 

Digits:=15;

Digits := 15

 (1)
> 

sola:=dsolve([diff(y(x), x$3)+diff(y(x), x$2)+y(x)=1, y(0)=0, y(1)=0, y(2)=1], [y(x)]):
> 

analoutput:=[evalf(subs(sola,x=0.5,y(x))),evalf(subs(sola,x=1.5,y(x)))];

analoutput := [-.126080336210879, .383074769960415]

 (2)
> 

eq1:=diff(y(x), x$3)+diff(y(x), x$2)+y(x)=1;

eq1 := diff(y(x), `$`(x, 3))+diff(y(x), `$`(x, 2))+y(x) = 1

 (3)
> 

eq2:=subs(y(x)=y2(x),eq1);

eq2 := diff(y2(x), `$`(x, 3))+diff(y2(x), `$`(x, 2))+y2(x) = 1

 (4)
> 

solnumeric:=dsolve([eq1,eq2,y(0)=0,y(1)=0,y2(0)=0,y2(1)=1,D(y)(1)=D(y2)(0),(D@@2)(y)(1)=(D@@2)(y2)(0)],[y(x),y2(x)],numeric):
> 

s1:=solnumeric(0.5);

s1 := [x = .5, y(x) = -.126080336419371, diff(y(x), x) = -0.433186250007679e-2, diff(y(x), `$`(x, 2)) = 1.01081632941025, y2(x) = .383074770066245, diff(y2(x), x) = 1.01399640036344, diff(y2(x), `$`(x, 2)) = .949493828849448]

 (5)
> 

numoutput:=[subs(s1,y(x)),subs(s1,y2(x))];

numoutput := [-.126080336419371, .383074770066245]

 (6)
> 

 


 

Download multipointBVP.mws

Dr. Venkat Subramanian

 

Try differentiating, there are no DAE so...

Dr. Venkat Subramanian 411
November 25 2016
2 4

Despite what Maple claims and what the help files say, there are no DAE solvers in Maple. Maple converts DAEs to dy/dt = f format and solves.

 

See my approach. See my past posts that helps avoid using fsolve. You have to make sure this code is correct( there may be more than one solution for initial condition).

In the future if possible please post y1, y2, etc instead of __t which took time for me to figure out.

> 

restart;
> 

Digits:=15;

Digits := 15

 (1)
> 

sys:={Q(0) = 0, Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)};

sys := {Q(0) = 0, Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, diff(Q(t), t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, diff(Q(t), t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t)}

 (2)
> 

indets(sys);

{t, (T__1(t)-T__1s(t))^1.32722881137112, (T__1(t)+T__1s(t))^.438263122318020, (T__2(t)+T__2s(t))^.438263121701134, (T__2s(t)-T__2(t))^1.32722881115416, Q(t), T__1(t), T__1s(t), T__2(t), T__2s(t), diff(Q(t), t)}

 (3)
> 

eqs:=[Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)];

eqs := [Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t), diff(Q(t), t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, diff(Q(t), t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (4)
> 

eqs[4]:=subs(eqs[3],eqs[4]);

eqs[4] := 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020

 (5)
> 

eqs[5]:=subs(eqs[3],eqs[5]);

eqs[5] := 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134

 (6)
> 

eqae:=[eqs[1],eqs[2],eqs[4],eqs[5]];

eqae := [Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (7)
> 

eqs1:=eval(subs(Q(t)=0,eqae));

eqs1 := [0 = 460900.000-5761.250*T__1(t), 0 = 5761.250*T__2(t)-220079.7500, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (8)
> 

eqs2:=subs(T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s,eqs1);

eqs2 := [0 = 460900.000-5761.250*T1, 0 = 5761.250*T2-220079.7500, 3300.00000000000*T1s-3300.00000000000*T2s = .654548222248727*(T1-T1s)^1.32722881137112*(T1+T1s)^.438263122318020, 3300.00000000000*T1s-3300.00000000000*T2s = .654548224676027*(T2s-T2)^1.32722881115416*(T2+T2s)^.438263121701134]

 (9)
> 

sol1:=fsolve({op(eqs2)});

sol1 := {T1 = 80.0000000000000, T1s = 60.3641019833216, T2 = 38.2000000000000, T2s = 60.2740143696713}

 (10)
> 

ics:=[Q(0)=0,T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s];

ics := [Q(0) = 0, T__1(t) = T1, T__2(t) = T2, T__1s(t) = T1s, T__2s(t) = T2s]

 (11)
> 

ics1:=subs(sol1,t=0,ics);

ics1 := [Q(0) = 0, T__1(0) = 80.0000000000000, T__2(0) = 38.2000000000000, T__1s(0) = 60.3641019833216, T__2s(0) = 60.2740143696713]

 (12)
> 

sys2:=diff(eqs[1],t),diff(eqs[2],t),diff(eqs[4],t),diff(eqs[5],t),eqs[3];

sys2 := diff(Q(t), t) = -5761.250*(diff(T__1(t), t)), diff(Q(t), t) = 5761.250*(diff(T__2(t), t)), 3300.00000000000*(diff(T__1s(t), t))-3300.00000000000*(diff(T__2s(t), t)) = .868735259000258*(T__1(t)-T__1s(t))^.32722881137112*(T__1(t)+T__1s(t))^.438263122318020*(diff(T__1(t), t)-(diff(T__1s(t), t)))+.286864347590436*(T__1(t)-T__1s(t))^1.32722881137112*(diff(T__1(t), t)+diff(T__1s(t), t))/(T__1(t)+T__1s(t))^.561736877681980, 3300.00000000000*(diff(T__1s(t), t))-3300.00000000000*(diff(T__2s(t), t)) = .868735262079829*(T__2s(t)-T__2(t))^.32722881115416*(T__2(t)+T__2s(t))^.438263121701134*(diff(T__2s(t), t)-(diff(T__2(t), t)))+.286864348250451*(T__2s(t)-T__2(t))^1.32722881115416*(diff(T__2(t), t)+diff(T__2s(t), t))/(T__2(t)+T__2s(t))^.561736878298866, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t)

 (13)
> 

 

> 

sol2:=dsolve({sys2,op(ics1)},type=numeric,implicit=true):
> 

sol2(1);

[t = 1., Q(t) = 296.806402013219, T__1(t) = 79.9484822908200, T__1s(t) = 60.3579593910879, T__2(t) = 38.2515177091800, T__2s(t) = 60.2681641420253]

 (14)
> 

plots:-odeplot(sol2,[t,Q(t)],0..1);

 
> 

plots:-odeplot(sol2,[t,T__1(t)],0..1);

 
> 

 


 

Download mapleprimesexample_DAE.mws
 

> 

restart;
> 

Digits:=15;

Digits := 15

 (1)
> 

sys:={Q(0) = 0, Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)};

sys := {Q(0) = 0, Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, diff(Q(t), t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, diff(Q(t), t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t)}

 (2)
> 

indets(sys);

{t, (T__1(t)-T__1s(t))^1.32722881137112, (T__1(t)+T__1s(t))^.438263122318020, (T__2(t)+T__2s(t))^.438263121701134, (T__2s(t)-T__2(t))^1.32722881115416, Q(t), T__1(t), T__1s(t), T__2(t), T__2s(t), diff(Q(t), t)}

 (3)
> 

eqs:=[Q(t) = (1.375*4190)*(80-T__1(t)), Q(t) = (1.375*4190)*(T__2(t)-38.2), diff(Q(t), t) = (240*0.1375e-1)*(T__1s(t)-T__2s(t))/(0.1e-2), diff(Q(t), t) = (0.1375e-1*47.6035070726347)*(T__1(t)-T__1s(t))*(T__1(t)+T__1s(t))^.438263122318020*((T__1(t)-T__1s(t))^.327228811371115), diff(Q(t), t) = (0.1375e-1*47.6035072491656)*(T__2s(t)-T__2(t))*(T__2(t)+T__2s(t))^.438263121701134*((T__2s(t)-T__2(t))^.327228811154163)];

eqs := [Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t), diff(Q(t), t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, diff(Q(t), t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (4)
> 

eqs[4]:=subs(eqs[3],eqs[4]);

eqs[4] := 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020

 (5)
> 

eqs[5]:=subs(eqs[3],eqs[5]);

eqs[5] := 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134

 (6)
> 

eqae:=[eqs[1],eqs[2],eqs[4],eqs[5]];

eqae := [Q(t) = 460900.000-5761.250*T__1(t), Q(t) = 5761.250*T__2(t)-220079.7500, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (7)
> 

eqs1:=eval(subs(Q(t)=0,eqae));

eqs1 := [0 = 460900.000-5761.250*T__1(t), 0 = 5761.250*T__2(t)-220079.7500, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548222248727*(T__1(t)-T__1s(t))^1.32722881137112*(T__1(t)+T__1s(t))^.438263122318020, 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t) = .654548224676027*(T__2s(t)-T__2(t))^1.32722881115416*(T__2(t)+T__2s(t))^.438263121701134]

 (8)
> 

eqs2:=subs(T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s,eqs1);

eqs2 := [0 = 460900.000-5761.250*T1, 0 = 5761.250*T2-220079.7500, 3300.00000000000*T1s-3300.00000000000*T2s = .654548222248727*(T1-T1s)^1.32722881137112*(T1+T1s)^.438263122318020, 3300.00000000000*T1s-3300.00000000000*T2s = .654548224676027*(T2s-T2)^1.32722881115416*(T2+T2s)^.438263121701134]

 (9)
> 

sol1:=fsolve({op(eqs2)});

sol1 := {T1 = 80.0000000000000, T1s = 60.3641019833216, T2 = 38.2000000000000, T2s = 60.2740143696713}

 (10)
> 

ics:=[Q(0)=0,T__1(t)=T1,T__2(t)=T2,T__1s(t)=T1s,T__2s(t)=T2s];

ics := [Q(0) = 0, T__1(t) = T1, T__2(t) = T2, T__1s(t) = T1s, T__2s(t) = T2s]

 (11)
> 

ics1:=subs(sol1,t=0,ics);

ics1 := [Q(0) = 0, T__1(0) = 80.0000000000000, T__2(0) = 38.2000000000000, T__1s(0) = 60.3641019833216, T__2s(0) = 60.2740143696713]

 (12)
> 

sys2:=diff(eqs[1],t),diff(eqs[2],t),diff(eqs[4],t),diff(eqs[5],t),eqs[3];

sys2 := diff(Q(t), t) = -5761.250*(diff(T__1(t), t)), diff(Q(t), t) = 5761.250*(diff(T__2(t), t)), 3300.00000000000*(diff(T__1s(t), t))-3300.00000000000*(diff(T__2s(t), t)) = .868735259000258*(T__1(t)-T__1s(t))^.32722881137112*(T__1(t)+T__1s(t))^.438263122318020*(diff(T__1(t), t)-(diff(T__1s(t), t)))+.286864347590436*(T__1(t)-T__1s(t))^1.32722881137112*(diff(T__1(t), t)+diff(T__1s(t), t))/(T__1(t)+T__1s(t))^.561736877681980, 3300.00000000000*(diff(T__1s(t), t))-3300.00000000000*(diff(T__2s(t), t)) = .868735262079829*(T__2s(t)-T__2(t))^.32722881115416*(T__2(t)+T__2s(t))^.438263121701134*(diff(T__2s(t), t)-(diff(T__2(t), t)))+.286864348250451*(T__2s(t)-T__2(t))^1.32722881115416*(diff(T__2(t), t)+diff(T__2s(t), t))/(T__2(t)+T__2s(t))^.561736878298866, diff(Q(t), t) = 3300.00000000000*T__1s(t)-3300.00000000000*T__2s(t)

 (13)
> 

 

> 

sol2:=dsolve({sys2,op(ics1)},type=numeric,implicit=true):
> 

sol2(1);

[t = 1., Q(t) = 296.806402013219, T__1(t) = 79.9484822908200, T__1s(t) = 60.3579593910879, T__2(t) = 38.2515177091800, T__2s(t) = 60.2681641420253]

 (14)
> 

plots:-odeplot(sol2,[t,Q(t)],0..1);

 
> 

plots:-odeplot(sol2,[t,T__1(t)],0..1);

 
> 

 


 

Download mapleprimesexample_DAE.mws
 

See code attached...

Dr. Venkat Subramanian 411
August 23 2016
0 0

p = df/deta, so, p, theta and F are solved.

 

Finite difference discretization is used.

 

Second order (central difference) for first derivative gives more oscillations. So first order approximation is used.

 

Code not very stable as fsolve does not work. Even NLPSolve obj = 1 is not robust. 

if you think these plots are meaningful. (sometimes P turns out to be negative), then you can build on this code for large values of inf (increae until there is no change) and increase N.

 

Check for typos and mistakes in the code as well.

> 

restart;
> 

Digits:=15;

Digits := 15

 (1)
> 

equ1 := (1+2*n)*f(eta)*(diff(theta(eta), eta))/(1+3*n)-(diff(theta(eta), `$`(eta, 2))) = 0;

equ1 := (1+2*n)*f(eta)*(diff(theta(eta), eta))/(1+3*n)-(diff(theta(eta), `$`(eta, 2))) = 0

 (2)
> 

equ2 := ((1+n)*(diff(f(eta), eta))^2/(1+3*n)-(1+2*n)*f(eta)*(diff(f(eta), eta, eta))/(1+3*n))/Bo+(diff(f(eta), `$`(eta, 3)))^n-theta(eta) = 0;

equ2 := ((1+n)*(diff(f(eta), eta))^2/(1+3*n)-(1+2*n)*f(eta)*(diff(f(eta), `$`(eta, 2)))/(1+3*n))/Bo+(diff(f(eta), `$`(eta, 3)))^n-theta(eta) = 0

 (3)
> 

Bo:=1;n:=2;

Bo := 1

n := 2

 (4)
> 

inf:=4;

inf := 4

 (5)
> 

bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(inf) = 0, theta(0) = 1, theta(inf) = 0;#(D@@3)(f)(0) = z(0);

bcs := f(0) = 0, (D(f))(0) = 0, (D(f))(4) = 0, theta(0) = 1, theta(4) = 0

 (6)
> 

equ1;equ2;

(5/7)*(f(eta)*(diff(theta(eta), eta)))-(diff(theta(eta), `$`(eta, 2))) = 0

(3/7)*((diff(f(eta), eta))^2)-(5/7)*(f(eta)*(diff(f(eta), `$`(eta, 2))))+(diff(f(eta), `$`(eta, 3)))^2-theta(eta) = 0

 (7)
> 

equ2:=subs(diff(f(eta),eta)=p(eta),equ2);

equ2 := (3/7)*(p(eta)^2)-(5/7)*(f(eta)*(diff(p(eta), eta)))+(diff(p(eta), `$`(eta, 2)))^2-theta(eta) = 0

 (8)
> 

equ3:=diff(f(eta),eta)-p(eta)=0;

equ3 := diff(f(eta), eta)-p(eta) = 0

 (9)
> 

N:=20;

N := 20

 (10)
> 

Eq1[0]:=T[0]=1;Eq2[0]:=P[0]=0;Eq3[0]:=F[0]=0;

Eq1[0] := T[0] = 1

Eq2[0] := P[0] = 0

Eq3[0] := F[0] = 0

 (11)
> 

#for i from 1 to N do Eq1[i]:=subs(diff(theta(eta),eta$2)=(T[i+1]-2*T[i]+T[i-1])/h^2,diff(theta(eta),eta)=(T[i+1]-T[i-1])/2/h,f(eta)=F[i],equ1);od;
> 

for i from 1 to N do Eq1[i]:=subs(diff(theta(eta),eta$2)=(T[i+1]-2*T[i]+T[i-1])/h^2,diff(theta(eta),eta)=(T[i]-T[i-1])/1/h,f(eta)=F[i],equ1);od;

Eq1[1] := (5/7)*(F[1]*(T[1]-T[0])/h)-(T[2]-2*T[1]+T[0])/h^2 = 0

Eq1[2] := (5/7)*(F[2]*(T[2]-T[1])/h)-(T[3]-2*T[2]+T[1])/h^2 = 0

Eq1[3] := (5/7)*(F[3]*(T[3]-T[2])/h)-(T[4]-2*T[3]+T[2])/h^2 = 0

Eq1[4] := (5/7)*(F[4]*(T[4]-T[3])/h)-(T[5]-2*T[4]+T[3])/h^2 = 0

Eq1[5] := (5/7)*(F[5]*(T[5]-T[4])/h)-(T[6]-2*T[5]+T[4])/h^2 = 0

Eq1[6] := (5/7)*(F[6]*(T[6]-T[5])/h)-(T[7]-2*T[6]+T[5])/h^2 = 0

Eq1[7] := (5/7)*(F[7]*(T[7]-T[6])/h)-(T[8]-2*T[7]+T[6])/h^2 = 0

Eq1[8] := (5/7)*(F[8]*(T[8]-T[7])/h)-(T[9]-2*T[8]+T[7])/h^2 = 0

Eq1[9] := (5/7)*(F[9]*(T[9]-T[8])/h)-(T[10]-2*T[9]+T[8])/h^2 = 0

Eq1[10] := (5/7)*(F[10]*(T[10]-T[9])/h)-(T[11]-2*T[10]+T[9])/h^2 = 0

Eq1[11] := (5/7)*(F[11]*(T[11]-T[10])/h)-(T[12]-2*T[11]+T[10])/h^2 = 0

Eq1[12] := (5/7)*(F[12]*(T[12]-T[11])/h)-(T[13]-2*T[12]+T[11])/h^2 = 0

Eq1[13] := (5/7)*(F[13]*(T[13]-T[12])/h)-(T[14]-2*T[13]+T[12])/h^2 = 0

Eq1[14] := (5/7)*(F[14]*(T[14]-T[13])/h)-(T[15]-2*T[14]+T[13])/h^2 = 0

Eq1[15] := (5/7)*(F[15]*(T[15]-T[14])/h)-(T[16]-2*T[15]+T[14])/h^2 = 0

Eq1[16] := (5/7)*(F[16]*(T[16]-T[15])/h)-(T[17]-2*T[16]+T[15])/h^2 = 0

Eq1[17] := (5/7)*(F[17]*(T[17]-T[16])/h)-(T[18]-2*T[17]+T[16])/h^2 = 0

Eq1[18] := (5/7)*(F[18]*(T[18]-T[17])/h)-(T[19]-2*T[18]+T[17])/h^2 = 0

Eq1[19] := (5/7)*(F[19]*(T[19]-T[18])/h)-(T[20]-2*T[19]+T[18])/h^2 = 0

Eq1[20] := (5/7)*(F[20]*(T[20]-T[19])/h)-(T[21]-2*T[20]+T[19])/h^2 = 0

 (12)
> 

for i from 1 to N do #Eq2[i]:=subs(diff(p(eta),eta$2)=(P[i+1]-2*P[i]+P[i-1])/h^2,diff(p(eta),eta)=(P[i+1]-P[i-1])/2/h,f(eta)=F[i],theta(eta)=T[i],p(eta)=P[i],equ2);od;
Eq2[i]:=subs(diff(p(eta),eta$2)=(P[i+1]-2*P[i]+P[i-1])/h^2,diff(p(eta),eta)=(P[i]-P[i-1])/1/h,f(eta)=F[i],theta(eta)=T[i],p(eta)=P[i],equ2);od;

Eq2[1] := (3/7)*(P[1]^2)-(5/7)*(F[1]*(P[1]-P[0])/h)+(P[2]-2*P[1]+P[0])^2/h^4-T[1] = 0

Eq2[2] := (3/7)*(P[2]^2)-(5/7)*(F[2]*(P[2]-P[1])/h)+(P[3]-2*P[2]+P[1])^2/h^4-T[2] = 0

Eq2[3] := (3/7)*(P[3]^2)-(5/7)*(F[3]*(P[3]-P[2])/h)+(P[4]-2*P[3]+P[2])^2/h^4-T[3] = 0

Eq2[4] := (3/7)*(P[4]^2)-(5/7)*(F[4]*(P[4]-P[3])/h)+(P[5]-2*P[4]+P[3])^2/h^4-T[4] = 0

Eq2[5] := (3/7)*(P[5]^2)-(5/7)*(F[5]*(P[5]-P[4])/h)+(P[6]-2*P[5]+P[4])^2/h^4-T[5] = 0

Eq2[6] := (3/7)*(P[6]^2)-(5/7)*(F[6]*(P[6]-P[5])/h)+(P[7]-2*P[6]+P[5])^2/h^4-T[6] = 0

Eq2[7] := (3/7)*(P[7]^2)-(5/7)*(F[7]*(P[7]-P[6])/h)+(P[8]-2*P[7]+P[6])^2/h^4-T[7] = 0

Eq2[8] := (3/7)*(P[8]^2)-(5/7)*(F[8]*(P[8]-P[7])/h)+(P[9]-2*P[8]+P[7])^2/h^4-T[8] = 0

Eq2[9] := (3/7)*(P[9]^2)-(5/7)*(F[9]*(P[9]-P[8])/h)+(P[10]-2*P[9]+P[8])^2/h^4-T[9] = 0

Eq2[10] := (3/7)*(P[10]^2)-(5/7)*(F[10]*(P[10]-P[9])/h)+(P[11]-2*P[10]+P[9])^2/h^4-T[10] = 0

Eq2[11] := (3/7)*(P[11]^2)-(5/7)*(F[11]*(P[11]-P[10])/h)+(P[12]-2*P[11]+P[10])^2/h^4-T[11] = 0

Eq2[12] := (3/7)*(P[12]^2)-(5/7)*(F[12]*(P[12]-P[11])/h)+(P[13]-2*P[12]+P[11])^2/h^4-T[12] = 0

Eq2[13] := (3/7)*(P[13]^2)-(5/7)*(F[13]*(P[13]-P[12])/h)+(P[14]-2*P[13]+P[12])^2/h^4-T[13] = 0

Eq2[14] := (3/7)*(P[14]^2)-(5/7)*(F[14]*(P[14]-P[13])/h)+(P[15]-2*P[14]+P[13])^2/h^4-T[14] = 0

Eq2[15] := (3/7)*(P[15]^2)-(5/7)*(F[15]*(P[15]-P[14])/h)+(P[16]-2*P[15]+P[14])^2/h^4-T[15] = 0

Eq2[16] := (3/7)*(P[16]^2)-(5/7)*(F[16]*(P[16]-P[15])/h)+(P[17]-2*P[16]+P[15])^2/h^4-T[16] = 0

Eq2[17] := (3/7)*(P[17]^2)-(5/7)*(F[17]*(P[17]-P[16])/h)+(P[18]-2*P[17]+P[16])^2/h^4-T[17] = 0

Eq2[18] := (3/7)*(P[18]^2)-(5/7)*(F[18]*(P[18]-P[17])/h)+(P[19]-2*P[18]+P[17])^2/h^4-T[18] = 0

Eq2[19] := (3/7)*(P[19]^2)-(5/7)*(F[19]*(P[19]-P[18])/h)+(P[20]-2*P[19]+P[18])^2/h^4-T[19] = 0

Eq2[20] := (3/7)*(P[20]^2)-(5/7)*(F[20]*(P[20]-P[19])/h)+(P[21]-2*P[20]+P[19])^2/h^4-T[20] = 0

 (13)
> 

#for i from 1 to N do Eq3[i]:=subs(diff(f(eta),eta)=(F[i+1]-F[i-1])/2/h,p(eta)=P[i],equ3);od;
> 

for i from 1 to N do Eq3[i]:=subs(diff(f(eta),eta)=(F[i]-F[i-1])/1/h,p(eta)=P[i],equ3);od;

Eq3[1] := (F[1]-F[0])/h-P[1] = 0

Eq3[2] := (F[2]-F[1])/h-P[2] = 0

Eq3[3] := (F[3]-F[2])/h-P[3] = 0

Eq3[4] := (F[4]-F[3])/h-P[4] = 0

Eq3[5] := (F[5]-F[4])/h-P[5] = 0

Eq3[6] := (F[6]-F[5])/h-P[6] = 0

Eq3[7] := (F[7]-F[6])/h-P[7] = 0

Eq3[8] := (F[8]-F[7])/h-P[8] = 0

Eq3[9] := (F[9]-F[8])/h-P[9] = 0

Eq3[10] := (F[10]-F[9])/h-P[10] = 0

Eq3[11] := (F[11]-F[10])/h-P[11] = 0

Eq3[12] := (F[12]-F[11])/h-P[12] = 0

Eq3[13] := (F[13]-F[12])/h-P[13] = 0

Eq3[14] := (F[14]-F[13])/h-P[14] = 0

Eq3[15] := (F[15]-F[14])/h-P[15] = 0

Eq3[16] := (F[16]-F[15])/h-P[16] = 0

Eq3[17] := (F[17]-F[16])/h-P[17] = 0

Eq3[18] := (F[18]-F[17])/h-P[18] = 0

Eq3[19] := (F[19]-F[18])/h-P[19] = 0

Eq3[20] := (F[20]-F[19])/h-P[20] = 0

 (14)
> 

Eq1[N+1]:=T[N+1]=0;Eq2[N+1]:=P[N+1]=0;
#Eq3[N+1]:=3*F[N+1]-4*F[N]+F[N-1]=0;
Eq3[N+1]:=F[N+1]-F[N]=0;

Eq1[21] := T[21] = 0

Eq2[21] := P[21] = 0

Eq3[21] := F[21]-F[20] = 0

 (15)
> 

eqs:=seq(Eq1[i],i=0..N+1),seq(Eq2[i],i=0..N+1),seq(Eq3[i],i=0..N+1),h=inf/(N+1):
> 

guess:=[seq(T[i]=0.9,i=0..N+1),seq(P[i]=-0.1,i=0..N+1),seq(F[i]=0,i=0..N+1),h=inf/(N+1)];

guess := [T[0] = .9, T[1] = .9, T[2] = .9, T[3] = .9, T[4] = .9, T[5] = .9, T[6] = .9, T[7] = .9, T[8] = .9, T[9] = .9, T[10] = .9, T[11] = .9, T[12] = .9, T[13] = .9, T[14] = .9, T[15] = .9, T[16] = .9, T[17] = .9, T[18] = .9, T[19] = .9, T[20] = .9, T[21] = .9, P[0] = -.1, P[1] = -.1, P[2] = -.1, P[3] = -.1, P[4] = -.1, P[5] = -.1, P[6] = -.1, P[7] = -.1, P[8] = -.1, P[9] = -.1, P[10] = -.1, P[11] = -.1, P[12] = -.1, P[13] = -.1, P[14] = -.1, P[15] = -.1, P[16] = -.1, P[17] = -.1, P[18] = -.1, P[19] = -.1, P[20] = -.1, P[21] = -.1, F[0] = 0, F[1] = 0, F[2] = 0, F[3] = 0, F[4] = 0, F[5] = 0, F[6] = 0, F[7] = 0, F[8] = 0, F[9] = 0, F[10] = 0, F[11] = 0, F[12] = 0, F[13] = 0, F[14] = 0, F[15] = 0, F[16] = 0, F[17] = 0, F[18] = 0, F[19] = 0, F[20] = 0, F[21] = 0, h = 4/21]

 (16)
> 

#sol:=fsolve({eqs},{op(guess)});#fsolve may not work
> 

 

> 

sol1:=Optimization:-NLPSolve(1,[eqs],initialpoint=guess);sol:=sol1[2];

sol1 := [1., [h = .190476190476190, F[0] = -0.245848511932597e-306, F[1] = 0.125531794137643e-1, F[2] = 0.444900932976673e-1, F[3] = .102568707625894, F[4] = .180089527713489, F[5] = .270496559574267, F[6] = .367423438879149, F[7] = .464741796622332, F[8] = .556611251552454, F[9] = .637530657145949, F[10] = .702390697389549, F[11] = .746528774293237, F[12] = .774101055141526, F[13] = .789402612662973, F[14] = .796849441131466, F[15] = .800957725636523, F[16] = .806320570644696, F[17] = .808293616187705, F[18] = .810716675579060, F[19] = .809850051322689, F[20] = .808350391184117, F[21] = .808350391184117, P[0] = -0.251136715241999e-306, P[1] = 0.659041919222627e-1, P[2] = .167668797890491, P[3] = .304912725223193, P[4] = .406984305459871, P[5] = .474636917269085, P[6] = .508866116350634, P[7] = .510921378151713, P[8] = .482314638383140, P[9] = .424826879365849, P[10] = .340515211278902, P[11] = .231724903744361, P[12] = .144754474453517, P[13] = 0.803331769875969e-1, P[14] = 0.390958494595872e-1, P[15] = 0.215684936515533e-1, P[16] = 0.281549362929053e-1, P[17] = 0.103584891007983e-1, P[18] = 0.127210618046135e-1, P[19] = -0.454977734594926e-2, P[20] = -0.787321572749885e-2, P[21] = 0.112206673565123e-306, T[0] = 1., T[1] = .975698861811284, T[2] = .951356219330371, T[3] = .926866229035161, T[4] = .902034482732556, T[5] = .876594309564110, T[6] = .850217880706270, T[7] = .822522905155056, T[8] = .793076771414036, T[9] = .761400698985263, T[10] = .726977079988965, T[11] = .689263824291512, T[12] = .647720088247062, T[13] = .601800971262948, T[14] = .550950062317141, T[15] = .494586157674124, T[16] = .432080061892214, T[17] = .362716829944366, T[18] = .285725590162289, T[19] = .200242094330574, T[20] = .105339712310215, T[21] = 0.]]

sol := [h = .190476190476190, F[0] = -0.245848511932597e-306, F[1] = 0.125531794137643e-1, F[2] = 0.444900932976673e-1, F[3] = .102568707625894, F[4] = .180089527713489, F[5] = .270496559574267, F[6] = .367423438879149, F[7] = .464741796622332, F[8] = .556611251552454, F[9] = .637530657145949, F[10] = .702390697389549, F[11] = .746528774293237, F[12] = .774101055141526, F[13] = .789402612662973, F[14] = .796849441131466, F[15] = .800957725636523, F[16] = .806320570644696, F[17] = .808293616187705, F[18] = .810716675579060, F[19] = .809850051322689, F[20] = .808350391184117, F[21] = .808350391184117, P[0] = -0.251136715241999e-306, P[1] = 0.659041919222627e-1, P[2] = .167668797890491, P[3] = .304912725223193, P[4] = .406984305459871, P[5] = .474636917269085, P[6] = .508866116350634, P[7] = .510921378151713, P[8] = .482314638383140, P[9] = .424826879365849, P[10] = .340515211278902, P[11] = .231724903744361, P[12] = .144754474453517, P[13] = 0.803331769875969e-1, P[14] = 0.390958494595872e-1, P[15] = 0.215684936515533e-1, P[16] = 0.281549362929053e-1, P[17] = 0.103584891007983e-1, P[18] = 0.127210618046135e-1, P[19] = -0.454977734594926e-2, P[20] = -0.787321572749885e-2, P[21] = 0.112206673565123e-306, T[0] = 1., T[1] = .975698861811284, T[2] = .951356219330371, T[3] = .926866229035161, T[4] = .902034482732556, T[5] = .876594309564110, T[6] = .850217880706270, T[7] = .822522905155056, T[8] = .793076771414036, T[9] = .761400698985263, T[10] = .726977079988965, T[11] = .689263824291512, T[12] = .647720088247062, T[13] = .601800971262948, T[14] = .550950062317141, T[15] = .494586157674124, T[16] = .432080061892214, T[17] = .362716829944366, T[18] = .285725590162289, T[19] = .200242094330574, T[20] = .105339712310215, T[21] = 0.]

 (17)
> 

plot(subs(sol,[seq([i*h,T[i]],i=0..N+1)]),axes=boxed);

 
> 

plot(subs(sol,[seq([i*h,F[i]],i=0..N+1)]),axes=boxed);

 
> 

plot(subs(sol,[seq([i*h,P[i]],i=0..N+1)]),axes=boxed);

 
> 

 

> 

 

> 

 

> 

 

> 

 

> 

 

> 

 

> 

 

 

 BVPFD.mws

FWIW...

Dr. Venkat Subramanian 411
August 21 2016
1 2

My experience with optmization in Maple for the last 10 years or so.

 

Help pages are poorly documented. So, you have to try many things to figure out what you want to do.

NLP problems can be done with different approaches in Maple.

As of now, more reliable answers are arrived for sqp approach (bounded and constrained problems) only if you state the problem in algebraic form (i.e, variables y1 to yn, bounds, etc in explicit form). Depsite what the help page claims, the matrix vector f(Y) form or using functional form f(u1,u2..) is not reliable. There seems to be some  bug in finding the numerical jacobian or hessian that will make Maple suggest that "initial guess meets first order optimality conditions, etc). method = nonlinearsimplex is reasonably robust and can be used only for unbounded and uncostrained problems.

Also, GAMS has access to many sophisticated optimizers. Maple's nonlinear constrained optimization is limited to regular sqp. So, you will be limtied by the size of the problem you can handle. This is because, Maple stores the entire Jacobian, and does not use sparse storage or solver approaches. It is easy to verify. Just take d^y/dx^2 = y^2 with bcs, x =0, y =1 and x = 1, y =0. Convert to finite difference form, and give 1 as objective and give bounds of 0..1 to for y. You will hit the size limit N  (number of node points) very quickly (perhaps N = 200 or 500) is the maximum size of optimization variables Maple can handle.

I have also used the Global optimization package (there is no gurantee for global optimization for any commerical package, despite the name) currently offered. The current package includes Genetic algorithm, but the previous global optimization offered by Pinto (during Maple 12,13) was more robust and had a good reduced gradient approach which was fast.

Mathematica offers access to KNITRO which is very robust. I think GAMS provides access to SNOPT and KNITRO, may be even IPOPT. Of course, people routinely do MEX files from MATLAB for Ipopt, knitro, etc.

I wish Maple provides access to KNITRO, SNOPT and IPOPT. I will be surprised if it happens considering the direction Maple has taken recently (focussing more on elementary stuff, compared to high-end users). That is a shame as Maple has the ability to do analytic Jacobian.

In our group, we run IPOPT from Maple by writing the C files for objective, gradients, constraints, jacobians and hessians. Proper documentation and more useful help files for externalcalling will help me and many to create a dll or package that can be added to Maple in the future.

 

 

check past post and explanation given by...

Dr. Venkat Subramanian 411
August 18 2016
0 0

http://www.mapleprimes.com/questions/125273-Dsolve-Events-How-To-Control-For-A-Sign-Change#comment125661

I have a solution...

Dr. Venkat Subramanian 411
January 04 2016
1 1

Bendesarts

See below and attached.

 

The model you have might be from a Hamiltonian. This means that all the solvers in Maple will fail. Read about geometric numerical integration. Basically we need algorithms that conserve energy. If you know which one is like displacement, and which one is like velocity, you may be able to come up with symplectic Euler like explicit schemes.

 

Fortunately, Gauss implicit Runge Kutta works for these type of problems (symplectic integrators). The second order Gauss method is the implicit mid point method.

 

That is for a system of ODEs dY/dt = F, Y1 = Y0+h*F(Ymid) where Ymid is given by Ymid = (Y0+Y1)/2. If F is nonlinear, then one needs to do a Newton Raphson iteration at every time step. This is a A stable (not L stable ) method good for some stiff problems, but does not have stiff decay property. Also not algebraically stable.

 

More about the property of this method later. The code below works. You can see that both dt = 0.1 and dt = 0.2 give same results from 0 to 400 for time.

Few ways to improve the code. Let me know if there is interest

 

1. Write the code in hardware floats (Speed up). This will mean that fsolve should be avoided and Newton Raphson should be done without using LinearSolve or any linear algebra package. Read my past code/post on Newton Raphson. Convert the same to compiled form for even more speed. If not for speed, do this for robustness as fsolve will fail and crash for large systems of equations.

2. Write adaptive time stepping instead of constant dt. This wil speed up the process and will also guarantee that results meet the tolerance specified.

3. Use higher order methods  - Higher order algorithms or simply using extrapolation methods. This might compromise on stability (in particular extrapolation).

 

4. Write schemes other than implicit mid point to march forward in time.

5. Write this and related (Adatptive time stepping with events, hardware floats/compiled form) for generic systems of ODEs and DAEs. The code is posted for ease of learning, not optimized for efficiency or performance.

 

Hope this is useful

 

> 

restart;
> 

r:=sqrt((x(t)/a)^2+(z(t)/b)^2);
eqx:=diff(x(t),t)=alpha*(1-r^2)*x(t)+w*a/b*z(t);
eqz:=diff(z(t),t)=beta*(1-r^2)*z(t)-w*b/a*x(t);
EqSys:=[eqx,eqz];

params := alpha=1, beta=1, a=0.4, b=0.2, w=1;

EqSys := eval([eqx,eqz], [params]);
xmax := 0.8; zmax := 0.4;
tmax := 400;
ic:=[x(0)=0.4, z(0)=0];

r := (x(t)^2/a^2+z(t)^2/b^2)^(1/2)

eqx := diff(x(t), t) = alpha*(1-x(t)^2/a^2-z(t)^2/b^2)*x(t)+w*a*z(t)/b

eqz := diff(z(t), t) = beta*(1-x(t)^2/a^2-z(t)^2/b^2)*z(t)-w*b*x(t)/a

EqSys := [diff(x(t), t) = alpha*(1-x(t)^2/a^2-z(t)^2/b^2)*x(t)+w*a*z(t)/b, diff(z(t), t) = beta*(1-x(t)^2/a^2-z(t)^2/b^2)*z(t)-w*b*x(t)/a]

params := alpha = 1, beta = 1, a = .4, b = .2, w = 1

EqSys := [diff(x(t), t) = (1-6.250000000*x(t)^2-25.00000000*z(t)^2)*x(t)+2.000000000*z(t), diff(z(t), t) = (1-6.250000000*x(t)^2-25.00000000*z(t)^2)*z(t)-.5000000000*x(t)]

xmax := .8

zmax := .4

tmax := 400

ic := [x(0) = .4, z(0) = 0]

 (1)
> 

DEplot(EqSys, [x(t),z(t)], t= 0..tmax, [ic],linecolor=black, thickness=1,x(t)=-xmax..xmax, z(t)=-zmax..zmax, scaling=constrained,arrows=none);

DEplot([diff(x(t), t) = (1-6.250000000*x(t)^2-25.00000000*z(t)^2)*x(t)+2.000000000*z(t), diff(z(t), t) = (1-6.250000000*x(t)^2-25.00000000*z(t)^2)*z(t)-.5000000000*x(t)], [x(t), z(t)], t = 0 .. 400, [[x(0) = .4, z(0) = 0]], linecolor = black, thickness = 1, x(t) = -.8 .. .8, z(t) = -.4 .. .4, scaling = constrained, arrows = none)

 (2)
> 

sol:=dsolve({op(EqSys),op(ic)},type=numeric,range=0..400):
> 

with(plots):
> 

odeplot(sol,[x(t),z(t)],t=0..400,thickness=3,axes=boxed);

 
> 

sol(400);

[t = 400., x(t) = -.210132490266510, z(t) = .170179637817480]

 (3)
> 

r:='r':
> 

r:=sqrt((x(t)/a)^2+(z(t)/b)^2);
eqx:=diff(x(t),t)=alpha*(1-r^2)*x(t)+w*a/b*z(t);

r := (x(t)^2/a^2+z(t)^2/b^2)^(1/2)

eqx := diff(x(t), t) = alpha*(1-x(t)^2/a^2-z(t)^2/b^2)*x(t)+w*a*z(t)/b

 (4)
> 

eqz:=diff(z(t),t)=beta*(1-r(t)^2)*z(t)-w*b/a*x(t);

eqz := diff(z(t), t) = beta*(1-(x(t))(t)^2/a(t)^2-(z(t))(t)^2/b(t)^2)*z(t)-w*b*x(t)/a

 (5)
> 

ic2:=x(0)=0.4,z(0)=0;

ic2 := x(0) = .4, z(0) = 0

 (6)
> 

parss:=[alpha=1, beta=1, a=0.4, b=0.2, w=1];
assign(parss);

parss := [alpha = 1, beta = 1, a = .4, b = .2, w = 1]

 (7)

 

 
> 

Digits:=15;

Digits := 15

 (8)
> 

eq1:=(lhs-rhs)(subs(diff(x(t),t)=(xnew-xold)/dt,diff(z(t),t)=(znew-zold)/dt,x(t)=(xnew+xold)/2,z(t)=(znew+zold)/2,eqx));

eq1 := (xnew-xold)/dt-(1-6.25000000000000*(xnew/2+xold/2)^2-25.0000000000000*(znew/2+zold/2)^2)*(xnew/2+xold/2)-1.00000000000000*znew-1.00000000000000*zold

 (9)
> 

eq2:=(lhs-rhs)(subs(diff(x(t),t)=(xnew-xold)/dt,diff(z(t),t)=(znew-zold)/dt,x(t)=(xnew+xold)/2,z(t)=(znew+zold)/2,eqz));

eq2 := (znew-zold)/dt-(1-6.25000000000000*(xnew/2+xold/2)(t)^2-25.0000000000000*(znew/2+zold/2)(t)^2)*(znew/2+zold/2)+.250000000000000*xnew+.250000000000000*xold

 (10)
> 

ss:=proc(Y0,dt)
local xold,zold,eq1,eq2,s2;
xold:=Y0[1];zold:=Y0[2];
eq1:=(xnew-xold)/dt-(1-6.250000000*(1/2*xnew+1/2*xold)^2-25.00000000*(1/2*znew+1/2*zold)^2)*(1/2*xnew+1/2*xold)-1.00000000000000*znew-1.00000000000000*zold;
eq2:=(znew-zold)/dt-(1-6.250000000*(1/2*xnew+1/2*xold)^2-25.00000000*(1/2*znew+1/2*zold)^2)*(1/2*znew+1/2*zold)+.250000000000000*xnew+.250000000000000*xold;
s2:=fsolve({eq1,eq2},{xnew=xold,znew=zold});
Vector(2,[subs(s2,xnew),subs(s2,znew)]);
end proc;

ss := proc (Y0, dt) local xold, zold, eq1, eq2, s2; xold := Y0[1]; zold := Y0[2]; eq1 := (xnew-xold)/dt-(1-6.250000000*((xnew/2+(1/2)*xold)^2)-25.00000000*((znew/2+(1/2)*zold)^2))*(xnew/2+(1/2)*xold)-1.00000000000000*znew-1.00000000000000*zold; eq2 := (znew-zold)/dt-(1-6.250000000*((xnew/2+(1/2)*xold)^2)-25.00000000*((znew/2+(1/2)*zold)^2))*(znew/2+(1/2)*zold)+.250000000000000*xnew+.250000000000000*xold; s2 := fsolve({eq1, eq2}, {xnew = xold, znew = zold}); Vector(2, [subs(s2, xnew), subs(s2, znew)]) end proc

 (11)
> 

Y[0]:=Vector(2,[0.4,0]);

Y[0] := Vector(2, {(1) = .4, (2) = 0})

 (12)
> 

ss(Y[0],0.1);

Vector(2, {(1) = .398095042153290, (2) = -0.199546389616883e-1})

 (13)
> 

tmax:=400.;N:=4000;dt:=tmax/N;

tmax := 400.

N := 4000

dt := .100000000000000

 (14)
> 

for i from 1 to N do Y[i]:=ss(Y[i-1],dt);od:
> 

p1:=plot([seq([Y[i][1],Y[i][2]],i=0..N)],axes=boxed):
> 

N:=2000;dt:=tmax/N;

N := 2000

dt := .200000000000000

 (15)
> 

 

> 

for i from 1 to N do Y[i]:=ss(Y[i-1],dt);od:
> 

p2:=plot([seq([Y[i][1],Y[i][2]],i=0..N)],axes=boxed,color=blue,thickness=3,axes=boxed):
> 

display({p1,p2});

 
> 

 

 

Download Gaussmethodproceduralcode.mws

False transient and robust transient met...

Dr. Venkat Subramanian 411
August 23 2015
2 2

Al86

See code attached below. y and z are used instead of x and y. phi is used instead of u. Basically artificial time derivatives are added and finite difference or other methods are used in x and y  in transient method of lines. For your problem, that may not work (if you have sqrt or stiff problems).

Attached code does both standard false transient and new a robust false transient method of lines we proposed. You can read more about it at http://depts.washington.edu/maple/falsetransient_new.html

 

Hope you find the approach useful. For higher number of node points, using stiff =true and compile =true will help. Stady state is reached in t = 0.05. I think the profile makes sense. You might want to increase the number of node points.

 

> 

restart;
> 

with(plots);

[animate, animate3d, animatecurve, arrow, changecoords, complexplot, complexplot3d, conformal, conformal3d, contourplot, contourplot3d, coordplot, coordplot3d, densityplot, display, dualaxisplot, fieldplot, fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d, inequal, interactive, interactiveparams, intersectplot, listcontplot, listcontplot3d, listdensityplot, listplot, listplot3d, loglogplot, logplot, matrixplot, multiple, odeplot, pareto, plotcompare, pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d, polyhedra_supported, polyhedraplot, rootlocus, semilogplot, setcolors, setoptions, setoptions3d, spacecurve, sparsematrixplot, surfdata, textplot, textplot3d, tubeplot]

 (1)

Height of the domain (y-coordinate)

> 

h:=1;

h := 1

 (2)

Length of the domain (z-coordinate)

> 

L:=1;

L := 1

 (3)
> 

 

Governing Equation (Laplace's Equation)

> 

Eq1:=diff(phi(y,z),y$2)+diff(phi(y,z),z$2)-sqrt(phi(y,z))-diff(phi(y,z),y)^2/phi(y,z)^(3/2);

Eq1 := diff(phi(y, z), `$`(y, 2))+diff(phi(y, z), `$`(z, 2))-phi(y, z)^(1/2)-(diff(phi(y, z), y))^2/phi(y, z)^(3/2)

 (4)

Boundary condition at y=0

> 

bcy1:=phi(y,z)-1;#diff(phi(y,z),y);

bcy1 := phi(y, z)-1

 (5)

Boundary condition at y=1

> 

bcy2:=diff(phi(y,z),y);

bcy2 := diff(phi(y, z), y)

 (6)

Boundary condition at z=0

> 

bcz1:=phi(y,z)-1;

bcz1 := phi(y, z)-1

 (7)

Boundary condition at z=1

> 

bcz2:=diff(phi(y,z),z);

bcz2 := diff(phi(y, z), z)

 (8)

Finite difference scheme to be used (2nd order-central difference for 2nd derivative)

> 

d2phidy2:=(phi[m+1,n]-2*phi[m,n]+phi[m-1,n])/dely^2;
d2phidz2:=(phi[m,n+1]-2*phi[m,n]+phi[m,n-1])/delz^2;

d2phidy2 := (phi[m+1, n]-2*phi[m, n]+phi[m-1, n])/dely^2

d2phidz2 := (phi[m, n+1]-2*phi[m, n]+phi[m, n-1])/delz^2

 (9)

Finite difference scheme to be used (2nd order-central difference for 1st derivative)

> 

dphidy:=(phi[m+1,n]-phi[m-1,n])/(2*dely);
> 

dphidz:=(phi[m,n+1]-phi[m,n-1])/(2*delz);

dphidy := (1/2)*((phi[m+1, n]-phi[m-1, n])/dely)

dphidz := (1/2)*((phi[m, n+1]-phi[m, n-1])/delz)

 (10)

3 Point forward and backward differences for 1st derivative (to be applied at boundary conditions)

> 

dphidyf:=(-phi[2,n]+4*phi[1,n]-3*phi[0,n])/(2*dely);
dphidyb:=(phi[Numy-1,n]-4*phi[Numy,n]+3*phi[Numy+1,n])/(2*delz);
dphidzf:=(-phi[m,2]+4*phi[m,1]-3*phi[m,0])/(2*delz);
dphidzb:=(phi[m,Numz-1]-4*phi[m,Numz]+3*phi[m,Numz+1])/(2*delz);

dphidyf := (1/2)*((-phi[2, n]+4*phi[1, n]-3*phi[0, n])/dely)

dphidyb := (1/2)*((phi[Numy-1, n]-4*phi[Numy, n]+3*phi[Numy+1, n])/delz)

dphidzf := (1/2)*((-phi[m, 2]+4*phi[m, 1]-3*phi[m, 0])/delz)

dphidzb := (1/2)*((phi[m, Numz-1]-4*phi[m, Numz]+3*phi[m, Numz+1])/delz)

 (11)

Number of interior node points in y

> 

Numy:=5;

Numy := 5

 (12)

Number of interior node points in z

> 

Numz:=5;

Numz := 5

 (13)
> 

 

> 

dely:=(h)/(Numy+1);

dely := 1/6

 (14)
> 

delz:=(L)/(Numz+1);

delz := 1/6

 (15)
> 

 

Develop finite difference equations for z boundary conditions for all y

> 

for j from 0 to Numy+1 do
Eq[j,0]:=subs(diff(phi(y,z),z)=dphidzf,phi(y,z)=phi[j,0],m=j,bcz1):
Eq[j,Numz+1]:=subs(diff(phi(y,z),z)=dphidzb,phi(y,z)=phi[j,Numz+1],m=j,bcz2):
od;

Eq[0, 0] := phi[0, 0]-1

Eq[0, 6] := 3*phi[0, 4]-12*phi[0, 5]+9*phi[0, 6]

Eq[1, 0] := phi[1, 0]-1

Eq[1, 6] := 3*phi[1, 4]-12*phi[1, 5]+9*phi[1, 6]

Eq[2, 0] := phi[2, 0]-1

Eq[2, 6] := 3*phi[2, 4]-12*phi[2, 5]+9*phi[2, 6]

Eq[3, 0] := phi[3, 0]-1

Eq[3, 6] := 3*phi[3, 4]-12*phi[3, 5]+9*phi[3, 6]

Eq[4, 0] := phi[4, 0]-1

Eq[4, 6] := 3*phi[4, 4]-12*phi[4, 5]+9*phi[4, 6]

Eq[5, 0] := phi[5, 0]-1

Eq[5, 6] := 3*phi[5, 4]-12*phi[5, 5]+9*phi[5, 6]

Eq[6, 0] := phi[6, 0]-1

Eq[6, 6] := 3*phi[6, 4]-12*phi[6, 5]+9*phi[6, 6]

 (16)

Develop finite difference equations for y boundary conditions for all z

> 

for j from 0 to Numz+1 do
Eq[0,j]:=subs(diff(phi(y,z),y)=dphidyf,phi(y,z)=phi[0,n],n=j,bcy1):
Eq[Numy+1,j]:=subs(diff(phi(y,z),y)=dphidyb,phi(y,z)=phi[Numy+1,n],n=j,bcy2):
od;

Eq[0, 0] := phi[0, 0]-1

Eq[6, 0] := 3*phi[4, 0]-12*phi[5, 0]+9*phi[6, 0]

Eq[0, 1] := phi[0, 1]-1

Eq[6, 1] := 3*phi[4, 1]-12*phi[5, 1]+9*phi[6, 1]

Eq[0, 2] := phi[0, 2]-1

Eq[6, 2] := 3*phi[4, 2]-12*phi[5, 2]+9*phi[6, 2]

Eq[0, 3] := phi[0, 3]-1

Eq[6, 3] := 3*phi[4, 3]-12*phi[5, 3]+9*phi[6, 3]

Eq[0, 4] := phi[0, 4]-1

Eq[6, 4] := 3*phi[4, 4]-12*phi[5, 4]+9*phi[6, 4]

Eq[0, 5] := phi[0, 5]-1

Eq[6, 5] := 3*phi[4, 5]-12*phi[5, 5]+9*phi[6, 5]

Eq[0, 6] := phi[0, 6]-1

Eq[6, 6] := 3*phi[4, 6]-12*phi[5, 6]+9*phi[6, 6]

 (17)
> 

#printlevel:=2;
> 

 

> 

 

Develop finite difference equations for all interior points using the governing equation

> 

for j from 1 to Numy do
for k1 from 1 to Numz do
Eq[j,k1]:=subs(diff(phi(y,z),y$2)=d2phidy2,diff(phi(y,z),z$2)=d2phidz2,diff(phi(y,z),y)=dphidy,diff(phi(y,z),z)=dphidz,phi(y,z)=phi[j,k1],n=k1,m=j,Eq1);
od;od;
> 

#printlevel:=1;

Collect all equations into a single list

> 

eqs:=[seq(seq(Eq[p,q],p=0..Numy+1),q=0..Numz+1)]:

Collect all variables into a single list

> 

vars:=[seq(seq(phi[i,j],i=0..Numz+1),j=0..Numy+1)]:
> 

 

Count number of equations

> 

n1:=nops(eqs);

n1 := 49

 (18)
> 

 

Convert all variables from the form phi[i,j] to YY[i]

> 

vars2:=[seq(vars[i]=YY[i](t),i=1..n1)];

vars2 := [phi[0, 0] = YY[1](t), phi[1, 0] = YY[2](t), phi[2, 0] = YY[3](t), phi[3, 0] = YY[4](t), phi[4, 0] = YY[5](t), phi[5, 0] = YY[6](t), phi[6, 0] = YY[7](t), phi[0, 1] = YY[8](t), phi[1, 1] = YY[9](t), phi[2, 1] = YY[10](t), phi[3, 1] = YY[11](t), phi[4, 1] = YY[12](t), phi[5, 1] = YY[13](t), phi[6, 1] = YY[14](t), phi[0, 2] = YY[15](t), phi[1, 2] = YY[16](t), phi[2, 2] = YY[17](t), phi[3, 2] = YY[18](t), phi[4, 2] = YY[19](t), phi[5, 2] = YY[20](t), phi[6, 2] = YY[21](t), phi[0, 3] = YY[22](t), phi[1, 3] = YY[23](t), phi[2, 3] = YY[24](t), phi[3, 3] = YY[25](t), phi[4, 3] = YY[26](t), phi[5, 3] = YY[27](t), phi[6, 3] = YY[28](t), phi[0, 4] = YY[29](t), phi[1, 4] = YY[30](t), phi[2, 4] = YY[31](t), phi[3, 4] = YY[32](t), phi[4, 4] = YY[33](t), phi[5, 4] = YY[34](t), phi[6, 4] = YY[35](t), phi[0, 5] = YY[36](t), phi[1, 5] = YY[37](t), phi[2, 5] = YY[38](t), phi[3, 5] = YY[39](t), phi[4, 5] = YY[40](t), phi[5, 5] = YY[41](t), phi[6, 5] = YY[42](t), phi[0, 6] = YY[43](t), phi[1, 6] = YY[44](t), phi[2, 6] = YY[45](t), phi[3, 6] = YY[46](t), phi[4, 6] = YY[47](t), phi[5, 6] = YY[48](t), phi[6, 6] = YY[49](t)]

 (19)

Perturbation parameter to be used

> 

mu:=1e-3:

Substitute new variables into the equations

> 

Eqs:=subs(vars2,eqs):

Eqs2 is the standard false transient formulation

> 

Eqs2:=seq(diff(YY[i](t),t)=Eqs[i],i=1..n1):

Eqs3 is the perturbation approach described in the paper

> 

Eqs3:=seq(mu*diff(Eqs[i],t)=-Eqs[i],i=1..n1):
> 

 

This is an initial guess for all values of phi to be used in the IVP solver

> 

ics2:=seq(YY[i](0)=1,i=1..n1):

Solver the standard false transient formulation with Maple's dsolve

> 

temp:=time[real]():sol2a:=dsolve({Eqs2,ics2},type=numeric,implicit=true,stiff=true,compile=true):time[real]()-temp;

4.239

 (20)

Solver the perturbation formulation with Maple's dsolve

> 

temp:=time[real]():sol3a:=dsolve({Eqs3,ics2},type=numeric,implicit=true,stiff=true,compile=true):time[real]()-temp;

5.045

 (21)
> 

 

> 

 

Plot the convergence of the standard false transient (red) and the perturbation approach (green). YY[1] is phi at y=0, z=0

> 

t11:=time[real]():p2:=odeplot(sol2a,[t,YY[25](t)],0..10):t11-time[real]();
t11:=time[real]():p3:=odeplot(sol3a,[t,YY[25](t)],0..10,color=green):t11-time[real]();

display(p2,p3);

 

 

 

Warning, cannot evaluate the solution further right of 1.1718747, probably a singularity

-.344

-.250

  

Calculate converged value from false transient approach (value at pseudo time=50)

> 

s2:=sol2a(50);

Error, (in sol2a) cannot evaluate the solution further right of 1.1718747, probably a singularity

  

Calculate converged value from perturbation approach (value at pseudo time=10)

> 

s3:=sol3a(10);

s3 := [t = 10., YY[1](t) = 1.00000000000000, YY[2](t) = 1.00000000000000, YY[3](t) = 1.00000000000000, YY[4](t) = 1.00000000000000, YY[5](t) = 1.00000000000000, YY[6](t) = 1.00000000000000, YY[7](t) = 1.00000000000000, YY[8](t) = 1.00000000000000, YY[9](t) = .961447296206741, YY[10](t) = .937843026786737, YY[11](t) = .922893858043597, YY[12](t) = .913630297499805, YY[13](t) = .908526302510143, YY[14](t) = .906824970846922, YY[15](t) = 1.00000000000000, YY[16](t) = .936207764943202, YY[17](t) = .894340729755916, YY[18](t) = .866952797282561, YY[19](t) = .849711238249887, YY[20](t) = .840140157843807, YY[21](t) = .836949797708447, YY[22](t) = 1.00000000000000, YY[23](t) = .919001245433058, YY[24](t) = .864046370332304, YY[25](t) = .827346201134004, YY[26](t) = .803956658334986, YY[27](t) = .790887040936321, YY[28](t) = .786530501803432, YY[29](t) = 1.00000000000000, YY[30](t) = .907624912769414, YY[31](t) = .843932739224344, YY[32](t) = .800894773972370, YY[33](t) = .773249720631034, YY[34](t) = .757732069206587, YY[35](t) = .752559518731771, YY[36](t) = 1.00000000000000, YY[37](t) = .901071504452547, YY[38](t) = .832356447116507, YY[39](t) = .785652162300325, YY[40](t) = .755526643672720, YY[41](t) = .738574197103932, YY[42](t) = .732923381581003, YY[43](t) = 1.00000000000000, YY[44](t) = .898887035013591, YY[45](t) = .828497683080560, YY[46](t) = .780571291742977, YY[47](t) = .749618951353282, YY[48](t) = .732188239736381, YY[49](t) = .726378002530746]

 (22)
> 

odeplot(sol3a,[t,YY[9](t)],0..0.05);

 
> 

 

Convert the solution (from s3, as a list), as an array

> 

for p from 0 to Numy+1 do
> 

for q from 0 to Numz+1 do
phi3[p,q]:=subs(vars2,s3,phi[p,q]);
yt:=p/(Numy+1);
zt:=q/(Numz+1);
od;
od:
> 

 

> 

 

> 

 

Plot the solution as found from the perturbation approach

> 

sol3:=matrix([seq([seq(phi3[p,q],p=0..Numy+1)],q=0..Numz+1)]):
> 

matrixplot(sol3,axes=boxed);

 
> 

 

> 

 

> 

 

> 

 

> 

 

> 

 

> 

 


Download mapleprimeslaplace.mws

Try DAE form - there are 3 solutions, al...

Dr. Venkat Subramanian 411
August 09 2015
1 7

See below. The problem is converted to index 1 DAE and solved. This problem might have multiple solutions (bifurcation) and might become index 2 DAE. I will explore more. Even at t= 0, there are 3 possible solutions.

> 

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;
odes2:=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});

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

odes2 := {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)}

 (1)
> 

odes3:=subs(diff(x1(t),t$2)=z1(t),diff(x2(t),t$2)=z2(t),diff(x1(t),t)=x1t(t),diff(x2(t),t)=x2t(t),{odes});

odes3 := {z1(t)*x2(t)+z2(t)*x1t(t)*sin(x1(t)^2)+cos(z2(t)*x2(t)) = sin(t), x1t(t)*z2(t)*sin(x1(t)*x2(t))+5*z1(t)*x2t(t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2 = exp(-x2(t)^2)}

 (2)
> 

odes4:=odes3 union {diff(x1(t),t)=x1t(t),diff(x1t(t),t)=z1(t),diff(x2(t),t)=x2t(t),diff(x2t(t),t)=z2(t)};

odes4 := {z1(t)*x2(t)+z2(t)*x1t(t)*sin(x1(t)^2)+cos(z2(t)*x2(t)) = sin(t), x1t(t)*z2(t)*sin(x1(t)*x2(t))+5*z1(t)*x2t(t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2 = exp(-x2(t)^2), diff(x1(t), t) = x1t(t), diff(x1t(t), t) = z1(t), diff(x2(t), t) = x2t(t), diff(x2t(t), t) = z2(t)}

 (3)
> 

Digits:=15;

Digits := 15

 (4)
> 

eqic:=eval(subs(x2(t)=2,x1(t)=1,x1t(t)=1,x2t(t)=2,z1(t)=z10,z2(t)=z20,t=0.,odes3));

eqic := {z20*sin(2)+10*z10*cos(1) = exp(-4), 2*z10+z20*sin(1)+cos(2*z20) = 0.}

 (5)
> 

z10:=solve(eqic[1],z10);

z10 := (1/10)*((-z20*sin(2)+exp(-4))/cos(1))

 (6)
> 

eqic[2];

(1/5)*((-z20*sin(2)+exp(-4))/cos(1))+z20*sin(1)+cos(2*z20) = 0.

 (7)
> 

zs:=[solve(eqic[2],z20)];

zs := [1.78622077114739, 1.07686043504888, -.627724205824561]

 (8)
> 

z10:=subs(z20=alpha,z10);

z10 := (1/10)*((-alpha*sin(2)+exp(-4))/cos(1))

 (9)
> 

sol:=dsolve({op(odes4),x2(0)=2,x1(0)=1,x1t(0)=1,x2t(0)=2,z1(0)=z10,z2(0)=alpha},type=numeric,compile=true,'parameters'=[alpha]):
> 

infolevel[all]:=0;

infolevel[all] := 0

 (10)
> 

sol('parameters'=[zs[1]]);

[alpha = 1.78622077114739]

 (11)
> 

sol(0.2);

[t = .2, x1(t) = 1.19406738431157, x1t(t) = .938736187324859, x2(t) = 2.43474508139573, x2t(t) = 2.34495747046976, z1(t) = -.388912586806758, z2(t) = 1.73594489563606]

 (12)
> 

sol('parameters'=[zs[2]]);

[alpha = 1.07686043504888]

 (13)
> 

sol(0.2);

[t = .2, x1(t) = 1.19658860073288, x1t(t) = .964367009670570, x2(t) = 2.41869398436184, x2t(t) = 2.16988682273856, z1(t) = -.262632675237659, z2(t) = .497906483056320]

 (14)
> 

sol('parameters'=[zs[3]]);

[alpha = -.627724205824561]

 (15)
> 

sol(0.2);

[t = .2, x1(t) = 1.20183378152547, x1t(t) = 1.01189227271276, x2(t) = 2.38886235498530, x2t(t) = 1.89879501190660, z1(t) = -.173626034799799, z2(t) = -.233760865624970]

 (16)
> 

with(plots):
> 

odeplot(sol,[t,z1(t)],0..3);

Warning, cannot evaluate the solution further right of .21144081, probably a singularity

  
> 

odeplot(sol,[t,z2(t)],0..3);

Warning, cannot evaluate the solution further right of .21144081, probably a singularity

  
> 

 


Download mapleprimesdaeapproach.mws

 

The procedure form works for this proble...

Dr. Venkat Subramanian 411
August 01 2015
0 0

method =lsode,

see below.

 

> 

restart:
> 

ff:= x1 -> evalf(Int(exp(t), t= 0 .. x1, method = _d01ajc) + 1);

ff := proc (x1) options operator, arrow; evalf(Int(exp(t), t = 0 .. x1, method = _d01ajc)+1) end proc

 (1)
> 

 

 

 

 

 
> 

p:=proc(N,x,Y,F)
global ff1;
local z1;
z1:=ff(x);
F[1]:=z1;
end proc;

p := proc (N, x, Y, F) local z1; global p1; z1 := ff(x); F[1] := z1 end proc

 (2)
> 

ics:=Array([1]);

ics := Matrix(1, 1, {(1, 1) = 1})

 (3)
> 

#infolevel[all]:=10;
> 

sol:=dsolve(numeric,number=1,procedure=p,initial=ics,procvars=[y(x)],start=0,method=lsode):
> 

sol(0.1);

[x = .1, y(x) = 1.10517091243516]

 (4)
> 

sol(1);

[x = 1., y(x) = 2.71828211984656]

 (5)
> 

 

> 

 

 

Download procedureapproach2.mws

the only way I know is to use known = ff...

Dr. Venkat Subramanian 411
August 01 2015
0 0

As I mentioned in the other post you are refering to.

The other approach is to use the procedure form of dsolve where I have shown to integrate a right hand side (rhs) which calls a simple function of x*f(x)^2. 

p1 can be your procedure where you do the complicated function and then you can call than in the main procedure p. The difficulty with the procedure form is that compile=true may not work (if p1 is not compiled), events/parametric form is a pain (almost impossible).

The attached code shows that you get the same results compared to direct numerical integration.

 

> 

restart:
> 

p1:=proc(x,f)
x*f^2;
end proc;

p1 := proc (x, f) x*f^2 end proc

 (1)

 

 

 

 
> 

p:=proc(N,x,Y,F)
global p1;
local z1;
z1:=p1(x,Y[1]);
F[1]:=z1;
end proc;

p := proc (N, x, Y, F) local z1; global p1; z1 := p1(x, Y[1]); F[1] := z1 end proc

 (2)
> 

ics:=Array([1]);

ics := Matrix(1, 1, {(1, 1) = 1})

 (3)
> 

#infolevel[all]:=10;
> 

sol:=dsolve(numeric,number=1,procedure=p,initial=ics,procvars=[y(x)],start=0):
> 

sol(1);

[x = 1., y(x) = 2.00000213434401]

 (4)
> 

eq2:=diff(y(x),x)=y(x)^2*x,y(0)=1;

eq2 := diff(y(x), x) = y(x)^2*x, y(0) = 1

 (5)
> 

sol2:=dsolve({eq2,y(0)=1},type=numeric):
> 

sol2(1);

[x = 1., y(x) = 2.00000213434401]

 (6)
> 

 

 

Download procedureapproach.mws

1 2 Page 1 of 2