Advanced Search

Venkat Subramanian

431 Reputation

13 Badges

15 years, 265 days
Contact 119368

MaplePrimes Activity


These are replies submitted by

different compilers...

Dr. Venkat Subramanian 431
August 11 2015
0 0

Preben,

My understanding is that Maple's Compiler:-compile and dsolve numeric use different compilers. It is possible to download and install mingw and use that for dsolve/numeric by linking to it from Maple (Thanks to Allan Wittkopf). If no one is sucessful here, I will dig my old codes and post. 

Also, this might have to do with windows 10 moving away from win32 to sandbox appraoches. Hope that is not the case and that it is an easy fix.

 

My 64bit maple compiles in Maple18/windows 7 configuration. However there is no gain in precision (You can use only Digits:=15. So, you don't lose anything staying with Maple 32 bit installations for compile=true approaches).

 

use events or stop_cond...

Dr. Venkat Subramanian 431
August 11 2015
0 0

@maple fan 

You have to know and state what causes the solution to stop(singular). For the 3rd root, you can  z2(t) close to zero to stop the solution (i.e, finding the t value for which the solution becomes singular). See below

 

> 

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],stop_cond=[z2(t)^2-1e-12]):
> 

infolevel[all]:=0;

infolevel[all] := 0

 (10)
> 

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

[alpha = 1.78622077114739]

 (11)
> 

sol(1000);

Error, (in sol) cannot evaluate the solution further right of .25108026, probably a singularity

  
> 

sol(0.25108026);

[t = .25108026, x1(t) = 1.24137113727403, x1t(t) = .907748693082798, x2(t) = 2.55688662743661, x2t(t) = 2.44179005569345, z1(t) = -1.18945323427187, z2(t) = 2.63365197729750]

 (12)
> 

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

[alpha = 1.07686043504888]

 (13)
> 

sol(1000);

Error, (in sol) cannot evaluate the solution further right of .21829382, probably a singularity

  
> 

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

[alpha = -.627724205824561]

 (14)
> 

sol(1000);

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

[t = .210891515596686, x1(t) = 1.21284210373439, x1t(t) = 1.00928747023560, x2(t) = 2.40953222982418, x2t(t) = 1.89717845937073, z1(t) = -.328141387462293, z2(t) = -0.999999999277117e-6]

 (15)
> 

with(plots):
> 

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

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

  
> 

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

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

  
> 

 

 

Download mapleprimesdaeapproach.mws

can you state your objective...

Dr. Venkat Subramanian 431
August 10 2015
0 0

@maple fan 

You are welcome. In my post above, I showed how to get the 3 possible solutions (using 3 different valeus of alpha (second derivative of x2 at t = 0)).

 

Are you trying to do bifuraction analysis? If so you have to use AUTO or someother specific software for that kind of analysis

tricks to speed up...

Dr. Venkat Subramanian 431
August 09 2015
0 0

@adel-00 

Are you sure you need Digits:=40;

Setting Digits:=15; and adding compile=true in the dsolve command (maxfun=0,compile=true) helps me plot the solution faster even with abserr=1e-21 and relerr=1e-9. You can even go for stricter tolerances (relerr=1e-13 and abserr=1e-25 or 1e-29, etc)

restart;
> Digits:=15:
> epsilon:=0.01:Delta1:=20:Delta2:=20: N1:=1000:
> dsys :={diff(x1(t),t)=-Delta2*x2(t)+y1(t)*epsilon, diff(x2(t),t)=Delta1*x1(t)+y2(t)*epsilon,diff(y1(t),t)=Delta2*y2(t)+x1(t)*z(t), diff(y2(t),t)=z(t)*x2(t)-y1(t)*Delta2, diff(z(t),t)=-4*x2(t)*y2(t)-4*y1(t)*x1(t)}:
> #infolevel[all]:=10;
> res:=dsolve(dsys union {x1(0)=0,x2(0)=0.2,y1(0)=0,y2(0)=0,z(0)=-1},numeric,output=listprocedure,maxfun=0,compile=true,abserr=1e-21,relerr=1e-9):
>
>
>
> plots[odeplot](res,[[t,y2(t)]],0..2000,axes=boxed,color=black,linestyle=1,tickmarks=[3, 4],thickness=2,view=[0..2000,-0.5..0.5]);

 

Try DAE form and use stop condition (eve...

Dr. Venkat Subramanian 431
August 09 2015
0 0

Preben, in the other post, I posted a DAE format of the problem. You can use stop condition to stop the procedure when z2 becomes zero. I am not sure if this is what you want. My guess is that time step for this problem keeps reducing (to 1e-15 or lower) as Maple tries to integrate past singularity.

 

When z2 = 0...

Dr. Venkat Subramanian 431
August 09 2015
0 0

@ the second derivative for x2 disappears from the original equations, one might have to differentiate the algebraic equation (twice) again to move forward (index 2 DAE)

See below...

Dr. Venkat Subramanian 431
August 02 2015
0 0

@Axel Vogt 

 

Higher values for x require longer time for integration and stricter tolerance (As the expected answer is 1e-18).

But dsolve approach works. 

 

> 

restart:
> 

Digits:=15;

Digits := 15

 (1)
> 

 

> 

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)

 (2)
> 

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

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

  

 

 
> 

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(2000);
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(2000); rhs(z1[2]) end if end proc

 (3)
> 

YFD2(45);

0.114346412039923e-17

 (4)
> 

 

 

Download knownfunctiondsolve5.mws

YFD is matching with YFD2...

Dr. Venkat Subramanian 431
August 02 2015
0 0

@Axel Vogt 

for different values of x in the code I posted above. If you just use sol1, then you have to call for different parameters (x). sol1 is a parametric dsolve. 

 

int(f(x,y)dy, y = x.. infinity) can be written as an ODE from 

dI/dy = f(x,y)

integrated from I(y=x) = 0 to I(y=infinity) to get the integral. 

integral...

Dr. Venkat Subramanian 431
August 02 2015
0 0

@Axel Vogt 

is with respect to xi right? 

Initial condition...

Dr. Venkat Subramanian 431
August 01 2015
0 0

@Axel Vogt 

 

I think sol2 is for the original ODE and needs the IC. I am using 0.153...as used by the OP.

Thanks...

Dr. Venkat Subramanian 431
August 01 2015
0 0

@Axel Vogt 

Bogo, also for this specific problem, I am not sure if the constants and ICs should be looked at. Both Y(x) and YFD(x) are very close (almost same).

improved range...

Dr. Venkat Subramanian 431
July 30 2015
0 3

@bogo 

Bogo, you are welcome. Thanks for the acknowledgement. You can call sol1(300) instead of sol1(100) to get sol2 for longer times. Also, for some reason, in Maple 14 classic worksheet it won't plot from 1 to 49 in the same graph, but the code wil plot if you split and plot. You can save those plots separately and display together.

See code below.

 

> 

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

y11:=subs(sol1(100),y(xi));

y11 := .153496865945855

 (11)
> 

y11-subs(sol1(2),y(xi));

.124847684227819

 (12)
> 

YFD(2);

0.995300744084111e-1

 (13)

 

 
> 

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(300);
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(300); rhs(z1[2]) end if end proc

 (14)
> 

YFD(1);YFD2(1);

.153496865945853

.153496865945855

 (15)

 

 

 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

 (16)
> 

sol2:=dsolve({de2,IC},type=numeric,known='YFD2',method=lsode,abserr=1e-4,maxfun=0):
> 

sol2(1);

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

 (17)

 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]

2.247

 (18)
> 

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

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

0.

 (19)
> 

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

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

18.283

 (20)
> 

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

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

0.

 (21)
> 

with(plots):
> 

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

 
> 

odeplot(sol2,[x,Y(x)],7..11);

 
> 

odeplot(sol2,[x,Y(x)],11..15);

 
> 

odeplot(sol2,[x,Y(x)],15..19);

 
> 

odeplot(sol2,[x,Y(x)],19..49);

 

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

Syntax...

Dr. Venkat Subramanian 431
July 29 2015
0 0

@bogo 

You are welcome.

I think it is just the syntax to enable dsolve/numeric codes to read. Not sure, perhpas Allan Wittkopf or other Maplesoft developers can explain the need for the same.

I have posted a more efficient approach. See below.

Thanks...

Dr. Venkat Subramanian 431
July 12 2015

@Preben Alsholm 

Preben, that is perfect. The parameters epsilon(=mu), tinit (initialization time), the constant q inside the tanh function are tunable parameters. My experience suggests that, 

epsilon = 1e-4 is good for abserr =1 e-6,

tinit = 0.1 is a good start for total time of integration of 1. Typically, we don't need more than tinit = 200 if we have to integrate till t = 10,000 for the original DAE system.

 

initstep = 0.1 or 0.01 is a good number (greater than epsilon).

q= 1000 is good, but for very long integrations, q = 100 may be sufficient.

See explanation below...

Dr. Venkat Subramanian 431
July 10 2015

@Markiyan Hirnyk 

you wrote "You wrote "if you take a linear system of DAEs. The mu approach will theoretically go to the expected steady state irrespective of the value of mu. (When only algebraic equations are solved). For example take Az=b as the algebraic equation". Can you kindly explain that in detail?"

First 10 11 12 13 14 15 16 Page 12 of 18