Venkat Subramanian

256 Reputation

13 Badges

9 years, 258 days

MaplePrimes Activity


These are answers submitted by

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

 

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
 

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

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.

 

 

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

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

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

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

 

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

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

Bogo, 

Please see improved code. I replaced the integral procedure with a numerical dsolve solution in compiled parametric form. You can see that the modified version is much faster. 

I know of at least one way to improve the speed even further, I will try that soon.

Hope you find this useful.


restart:

Digits:=15;

Digits := 15

(1)

de := diff(Y(x), x)+10^5*(Y(x)^2-YFD(x)^2)/x^2 = 0;

de := diff(Y(x), x)+100000*(Y(x)^2-YFD(x)^2)/x^2 = 0

(2)

YFD:=proc (x)  local xi,g; g:=2;
if not type(evalf(x),'numeric') then
      'procname'(x)
   else Re(evalf((1/2)*g*(int((xi^2-x^2)^(1/2)*xi/(exp(xi)+1), xi = x .. infinity))/Pi^2));
end if;
end proc;

YFD := proc (x) local xi, g; g := 2; if not type(evalf(x), 'numeric') then ('procname')(x) else Re(evalf((1/2)*((g/Pi^2)*(int((xi^2-x^2)^(1/2)*xi/(exp(xi)+1), xi = x .. infinity))))) end if end proc

(3)

 

 

IC:=Y(1)=YFD(1);

IC := Y(1) = .153496865945853

(4)

#infolevel[all]:=10;

sol:=dsolve({de,IC},type=numeric,known='YFD',method=lsode):

sol(1);

[x = 1., Y(x) = .153496865945853]

(5)

 Solution will be evaluated at a later stage

#sol(1.00001);

#sol(1.0001);

 

 

 The itnegral equation is converted to a differential equation and solved in parametric and compiled form with respect to x and xi = 100 is taken to be infinity (not changing much beyond x = 45).

eq2:=diff(y(xi),xi)=(xi^2-1.*x^2)^(1/2)*xi/(exp(xi)+1.)/Pi^2;

eq2 := diff(y(xi), xi) = (xi^2-1.*x^2)^(1/2)*xi/((exp(xi)+1.)*Pi^2)

(6)

sol1:=dsolve({eq2,y(x)=0},type=numeric,'parameters'=[x],compile=true,abserr=1e-15,relerr=1e-14):

Warning, relerr increased from .1e-13 to 5.00001e-013

 

sol1('parameters'=[1]);

[x = 1.]

(7)

sol1(1);sol1(2);

[xi = 1., y(xi) = 0.]

[xi = 2., y(xi) = 0.286491817180358e-1]

(8)

sol1(45.00000000000);

[xi = 45.00000000000, y(xi) = .153496865945855]

(9)

sol1(100);

[xi = 100., y(xi) = .153496865945855]

(10)

 

 

YFD2:=proc (x)  
global sol1,y;
local xi,z1,g; g:=2;
if not type(evalf(x),'numeric') then
      'procname'(x)
   else
   sol1('parameters'=[x]);
z1:=sol1(100);
rhs(z1[2]);#subs(sol1(100),y(xi)):   

end if;
end proc;

YFD2 := proc (x) local xi, z1, g; global sol1, y; g := 2; if not type(evalf(x), 'numeric') then ('procname')(x) else sol1('parameters' = [x]); z1 := sol1(100); rhs(z1[2]) end if end proc

(11)

YFD(1);YFD2(1);

.153496865945853

.153496865945855

(12)

 

 

 One can note that the dsolve solution for the integral is much faster than the integral calucation and provides the same value.

de2:=diff(Y(x), x)+10^5*(Y(x)^2-YFD2(x)^2)/x^2 = 0;;

de2 := diff(Y(x), x)+100000*(Y(x)^2-YFD2(x)^2)/x^2 = 0

(13)

sol2:=dsolve({de2,IC},type=numeric,known='YFD2',method=lsode):

sol2(1);

[x = 1., Y(x) = .153496865945853]

(14)

 Next time taken is compared for the original dsolve with integral for the procedure and the new dsolve.

t11:=time():sol(1.00001);time()-t11;

[x = 1.00001, Y(x) = .153496759729280]

4.758

(15)

t11:=time():sol2(1.00001);time()-t11;

[x = 1.00001, Y(x) = .153496759729280]

0.

(16)

t11:=time():sol(1.0001);time()-t11;

[x = 1.0001, Y(x) = .153493467077101]

28.642

(17)

t11:=time():sol2(1.0001);time()-t11;

[x = 1.0001, Y(x) = .153493467077101]

0.46e-1

(18)

with(plots):

odeplot(sol2,[x,Y(x)],1..2);

 

 

 Note that in sol2, compile=true can't be used as it is not able to recognize YDF2. If that can be enabled this code will run even faster.


Download knownfunctiondsolve2.mws

see attached.

 

Read about known function for dsolve/numeric

http://www.maplesoft.com/support/help/Maple/view.aspx?path=dsolve/numeric

 

I used method lsode. As this problem is stiff, use stiff = true, but you have to provide the diff of the function (see the help above above).

 

restart:

de := diff(Y(x), x)+10^5*(Y(x)^2-YFD(x)^2)/x^2 = 0;

de := diff(Y(x), x)+100000*(Y(x)^2-YFD(x)^2)/x^2 = 0

(1)

YFD:=proc (x)  local xi; g:=2;
if not type(evalf(x),'numeric') then
      'procname'(x)
   else Re(evalf((1/2)*g*(int((xi^2-x^2)^(1/2)*xi/(exp(xi)+1), xi = x .. infinity))/Pi^2));
end if;
end proc;

Warning, `g` is implicitly declared local to procedure `YFD`

YFD := proc (x) local xi, g; g := 2; if not type(evalf(x), 'numeric') then ('procname')(x) else Re(evalf((1/2)*((g/Pi^2)*(int((xi^2-x^2)^(1/2)*xi/(exp(xi)+1), xi = x .. infinity))))) end if end proc

(2)

IC:=Y(1)=YFD(1);

IC := Y(1) = .1534968659

(3)

#infolevel[all]:=10;

sol:=dsolve({de,IC},type=numeric,known='YFD',method=lsode,maxfun=0):

sol(1);

[x = 1., Y(x) = .1534968659]

(4)

sol(1.00001);

[x = 1.00001, Y(x) = .153496759694509]

(5)

sol(1.0001);

[x = 1.0001, Y(x) = .153493467074749]

(6)

#sol(1.001);

 

 

Download knownfunctiondsolve.mws

 

I ran for shorttimes, it is taking long. You can read about lsode options if you want to use this. Also, please consider using classic worksheet or maple format for display. Most of the time was spent on copying and pasting from your code to mine, not the actual solution itself. Try classic worksheet, and see for yourself how easy it is to copy paste commands or write a code with large number of sentences/commands.

Without BCs.  May be you can see that and find out which one makes sense. A second order BVP can have only 2 BCs. you have 3 for each variable.

with(inttrans);

 

this should work

If you know what causes singularity, you can do this by finding the time at which event happens. See for example the code below. You can ofcourse use v(a) =b and call both initial t values and intial conditions as parameters in dsolve and then call the dsolve solution for different parameters.

 

 

restart;

Digits:=15;

Digits := 15

(1)

ode1:=
     -0.1*diff(u(z),z,z) +
     (z-2*diff(v(z)^(-1/2),z))*diff(u(z),z)+(3-2*diff(v(z)^(-1/2),z,z))*u(z)
          =
     0;

ode1 := -.1*(diff(u(z), `$`(z, 2)))+(z+(diff(v(z), z))/v(z)^(3/2))*(diff(u(z), z))+(3-(3/2)*((diff(v(z), z))^2/v(z)^(5/2))+(diff(v(z), `$`(z, 2)))/v(z)^(3/2))*u(z) = 0

(2)

ode2:= 0.1*diff(v(z),z,z)+0.01*z*diff(v(z),z)+0.02*v(z)-u(z)*v(z)^(1/2)=0;

ode2 := .1*(diff(v(z), `$`(z, 2)))+0.1e-1*z*(diff(v(z), z))+0.2e-1*v(z)-u(z)*v(z)^(1/2) = 0

(3)

ics:= u(0)=0, v(0)=0.1, D(u)(0)=0, D(v)(0)=0:

sol:= dsolve({ode||(1..2), ics}, numeric,stop_cond=[v(z)-1e-7],abserr=1e-10):

s1:=sol(100);

Warning, cannot evaluate the solution further right of 4.1328703, event #1 triggered a halt

s1 := [z = 4.13287033838952, u(z) = 0., diff(u(z), z) = 0., v(z) = 0.100000000000507e-6, diff(v(z), z) = -0.241962748394962e-1]

(4)

z1:=subs(s1,z);

z1 := 4.13287033838952

(5)

 


Download sinularity.mws

1 2 Page 1 of 2