Venkat Subramanian

386 Reputation

13 Badges

14 years, 111 days

MaplePrimes Activity


These are answers submitted by

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

If you know the specific integral you want to calculate beforehand, then add an additional variable z(t) defined as

dz/dt = y^2 with z(0)=0. z(1) gives the integral. The advantage is that the integral you get is of the same order of accuracy as u. You can use arbitrary number of Digits (Digits = 20 or 30).

Unfortunately this approach can be used only if the resulting system is an ODE. If your resulting system is a DAE, Maple cannot solve it (DAE BVP)

restart:

eq:=-1.2*diff(y(t),t$2)+0.2*y(t)=2;

eq := -1.2*(diff(y(t), `$`(t, 2)))+.2*y(t) = 2

(1)

bcs:=y(0)=1,y(1)=0;

bcs := y(0) = 1, y(1) = 0

(2)

eq2:=diff(z(t),t)=y(t)^2;

eq2 := diff(z(t), t) = y(t)^2

(3)

sol:=dsolve({eq,eq2,bcs,z(0)=0},type=numeric,abserr=1e-8,maxmesh=100):

with(plots):

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

 

sol(1);

[t = 1., y(t) = 0., diff(y(t), t) = -1.79470197426050, z(t) = .482858946230082]

(4)

 

 

Download IntegralMapleprimes.mws

hope this helps.

restart:

eq1:=diff(V(t)*C(t),t)=G-K*C(t);

eq1 := (diff(V(t), t))*C(t)+V(t)*(diff(C(t), t)) = G-K*C(t)

(1)

eq2:=diff(V(t),t)=alpha-beta;

eq2 := diff(V(t), t) = alpha-beta

(2)

dsolve({eq1,eq2});

[{V(t) = (alpha-beta)*t+_C2}, {C(t) = (Int(G*exp((alpha-beta+K)*(Int(1/V(t), t)))/V(t), t)+_C1)*exp(Int(-(alpha-beta+K)/V(t), t))}]

(3)

V(t):=c1+(alpha-beta)*t;

V(t) := c1+(alpha-beta)*t

(4)

eq1;

(alpha-beta)*C(t)+(c1+(alpha-beta)*t)*(diff(C(t), t)) = G-K*C(t)

(5)

dsolve({eq1,C(0)=c2});

C(t) = ((c1+t*alpha-t*beta)^(1+K/(alpha-beta))*G/(alpha-beta+K)-(-c2*alpha+c2*beta-c2*K+c1^(-(alpha-beta+K)/(alpha-beta))*c1^((alpha-beta+K)/(alpha-beta))*G)/(c1^(-(alpha-beta+K)/(alpha-beta))*(alpha-beta+K)))*(c1+(alpha-beta)*t)^(-alpha/(alpha-beta)+beta/(alpha-beta)-K/(alpha-beta))

(6)

 


Download solvesequentially.mws

Add u as an additional variable defined by dh/dt = u

Then enter to Maple in DAE form. My guess is that you are lucky because this approach works for this case as the addition of u enables Maple to convert this system to explicit dy/dt = f form.  But in general Maple cannot solve nonlinear DAEs or ODEs if it cannot convert them to explicit ODE form dy/dt = f. In fact, this is a serious flaw. The assumption made here is that all ODEs and DAEs can be converted to dy/dt = f(y) form. This is not true for many nonlinear or implicit ODEs until you add addtional variables. Maple avoids using Newton Raphson or any nonlinear solvers in its ODE solver.

Code attached. It is best to post some range for parameters or the physics of the problem. I can't guess the values of the parameters. You have to find u(0) for every set of parameters. For the parameters I chose, there is a singularity at 3.22 or so.

 

 

 

 

restart:

ODE1:=A*u(t)^2+B*u(t)*(h(t)+C)+E*h(t)=F*cos(u(t));

ODE1 := A*u(t)^2+B*u(t)*(h(t)+C)+E*h(t) = F*cos(u(t))

(1)

eqic:=A*u1^2+B*u1*C=F*cos(u1);

eqic := A*u1^2+B*C*u1 = F*cos(u1)

(2)

 For the simulation shown u has an initial condition of Pi/2. When parameters change, change ic accordingly.

ODE2:=diff(h(t),t)=u(t):

sol:=dsolve({ODE1,ODE2,h(0)=0,u(0)=Pi/2},type=numeric,parameters=[A,B,C,E,F],stiff=true,maxfun=0):

sol('parameters'=[0,0,1,1,1]):

with(plots):

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

 

odeplot(sol,[t,u(t)],0..1);

 

 

 

 

Download convertimplicitODEtoDAE.mws

I would suggest to post the original problem first so that we can make sure there was no mistake in the formulation when converted to adjoint formulation as a BVP.

If the code is posted as a mws for the BVP, I can try to solve and post it back. One can sometimes trick the BVPsolver in Maple to solve when the Newton iteration fails. Read help files on continuation.

 

1 2 Page 2 of 2