Venkat Subramanian

386 Reputation

13 Badges

15 years, 248 days

MaplePrimes Activity


These are replies submitted by

@JonMcLoone 

As mentioned before by me Maple does not have a direct DAE solver (it converts DAEs to ODEs and integrates as of Maple 2015). Mathematica has IDA (Backward difference based solver) for index 1 nonlinear DAEs and can solve DAEs Maple cannot solve. 

That was the only reason I tried to use Mathematica as I had to solve DAEs.

 

That said, my attempts to use Mathematica kept pushing me back to Maple (classic worksheet) as I couldn't handle 100s of lines (copy paste multiple procedures) in Mathematica (too slow). I would rank the regular window (mw) comparable to Mathematica in terms of speed and ease of use (copy/paste).

 

In addition, setting higher infolevel helps me see the algorithm behind a particular method or command in Maple (at least for most of it, not true for financial portfolio/newer packages). I don't think one can see the algorithm behind Mathematica. Please correct me if I am wrong. (For me this is very important. In fact I have learned more mathematial methods/algorithms from just browsing through Maple than reading text/papers. So, I have lost interest in Mathematica).

 

Mathematica used to provide access to KNITRO (not sure if it is true anymore) for optimization which got my attention for a while (KNITRO, SNOPT and IPOPT are some of the best approaches for solving large scale optimization problems. TOMLAB is a good start for those who want to see what is out there). The internal and external optimization package called by Maple is not the best and you will hit RAM limitations very fast (my guess is that zeros in Jacobian are stored and that consumes memory, moving to sparse storage would help).

 

When I teach, I am being pushed to use MATLAB by many, but I insist that students learn and program in what they like and I will use Maple for the symbolic stuff (generate equations/models) and use C/FORTRAN for performance.  MATLAB for control/optimization tool box.

 

 

 

@maple2015 

try  using fsolve first. If not, the past link might be useful.

http://www.mapleprimes.com/questions/204206-Towards-A-Compiled-Newton-Raphson

 

 

@Preben Alsholm 

Yes, theoretically fsolve should work. But efficiency of fsolve is pretty bad compared to what I posted. You can check the cpu time. Importantly for difficult problems and large scale problems fsolve may not work.

Also, one can specify abserr=1e-15, relerr=1e-13 to get more accurate solution for the dsolve solution. 

 

in a PDE form with well defined bcs and ics. Population balance models can be very generic and complicated (with moving particle size as opposed to constant particle size), x,y,z and time variables. So, you need to state the equations and then you can expect a Maple solution.

 

Of course there are some very simple models with hyperbolic/parabolic type equations with analytical solutions.

As far as I know, no one claim that they have a code to solve any population balance model with all the complexities (this includes FLUENT/CFD commerical software).

If you want to solve this, then search and try COLDAE (FORTRAN code).

For some problems you can differentiate the algebraic equations and use the equation as bc to solve this. If you want me to try this, post mws code or mws format code (Maple format as opposed math format). Your problem looks simple, so this approach might work

 

@tomleslie 

I checked both maple 18 and maple 2015 (not 2015.1) in 64 bit windows 7, dsolve numeric compile =true works. The 32 bit calls the watcom, but the 64 bit calls something else, I can post a screenshot for both 32 and 64 bit versions if that helps.

Regarding your point on the same compiler version, point well taken. It is trivial to write a compiled procedure (single one, not linking multiple procedures), but sadly dsolve numeric/compile won't recognize this.

 

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

 

@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

@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

@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]);

 

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.

 

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

@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

@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. 

@Axel Vogt 

is with respect to xi right? 

First 9 10 11 12 13 14 15 Page 11 of 17