Venkat Subramanian

386 Reputation

13 Badges

15 years, 331 days

MaplePrimes Activity


These are replies submitted by

@Markiyan Hirnyk 

@Markiyan Hirnyk 

At t= 0  the model is index 1 DAE. That is if you differentiate the second equation you get

 

2ydy/dt+2zdz/dt=0 

substituting first equation one gets

2zy+2zdz/dt=0

which implies z = 0 or dz/dt = -y

when z is not equal to zero one can solve 2 ODEs

dy/dt = 0 and dz/dt = - y

 

But when z = 0 we have

dy/dt = 0 and 

y^2- 1 = 0

 

which has two solutions for y, y =1 and y =-1. We need to differentiate the second equation twice to get an equation for dz/dt (hence this is called index 2 DAE).

 

 

When this DAE is sovled in standard index 1 DAE solvers in FORTRAN or MATLAB, the code stops at t = Pi/2 accurately. Maple's dsolve numeric will integrate past this singular point. Code attached for Maple's dsolve numeric. (Uncomment infolevel statement to see the equations being integrated by dsolve/numeric).

 

restart:

eq1:=diff(y(t),t)=z(t);

eq1 := diff(y(t), t) = z(t)

(1)

eq2:=y(t)^2+z(t)^2-1;

eq2 := y(t)^2+z(t)^2-1

(2)

#infolevel[all]:=10;

sol:=dsolve({eq1,eq2,y(0)=0,z(0)=1},type=numeric):

with(plots):

odeplot(sol,[[t,y(t)],[t,z(t)]],0..10);

 

 

 

 

Download index2example.mws

@Markiyan Hirnyk 

 

Markiyan, I use my middle initial as well in publications. My google schoalr page is https://scholar.google.com/citations?user=_rQBjg4AAAAJ&hl=en. The author you are refering to is some one else.

Thanks

 

Note that when dsolve starts, the time till t = tinit is for initialization. Setting initstep=0.1 (to avoid Maple spending too much time near t =0) and compile =true will help solve many problems efficiently (can't be done for Digits greater than 15 as of Maple 2015).

Modified version of example 4 is attached to show the improvement in the efficiency.

 

Example 4

restart;

with(plots):

N:=11:h:=1/(N+1):

for i from 2 to N+1 do eq1[i]:=diff(y||i(t),t)=(y||(i+1)(t)-2*y||i(t)+y||(i-1)(t))/h^2-y||i(t)*(1+z||i(t));od:

for i from 2 to N+1 do eq2[i]:=0=(z||(i+1)(t)-2*z||i(t)+z||(i-1)(t))/h^2-(1-y||i(t)^2)*(exp(-z||i(t)));od:

eq1[1]:=(3*y1(t)-4*y2(t)+y3(t))/(2*h)=0: eq1[N+2]:=y||(N+2)(t)-1=0:

eq2[1]:=(3*z1(t)-4*z2(t)+1*z3(t))/(2*h)=0: eq2[N+2]:=z||(N+2)(t)=0:

eq1[1]:=subs(y1(t)=z||(N+3)(t),eq1[1]):

eq1[N+2]:=subs(y||(N+2)(t)=z||(N+4)(t),eq1[N+2]):

eqode:=[seq(subs(y1(t)=z||(N+3)(t),y||(N+2)(t)=z||(N+4)(t),eq1[i]),i=2..N+1)]:

eqae:=[eq1[1],eq1[N+2],seq(eq2[i],i=1..N+2)]:

icodes:=[seq(y||j(0)=1,j=2..N+1)]:

icaes:=[seq(z||j(0)=0,j=1..N+2),z||(N+3)(0)=1,z||(N+4)(0)=1]:

pars:=[mu=0.00001,q=1000,tint=1,tf=2]:

Xexplicit:=2:

ff:=subs(pars,1/2+1/2*tanh(q*(t-tint))):

NODE:=nops(eqode):NAE:=nops(eqae):

for XX from 1 to NODE do

EQODE||XX:=lhs(eqode[XX])=rhs(eqode[XX])*ff: end do:

for XX from 1 to NAE do

EQAE||XX:=subs(pars,-mu*(diff(rhs(eqae[XX])-lhs(eqae[XX]),t))=rhs(eqae[XX])-lhs(eqae[XX])); end do:

Dvars1:={seq(diff(z||x(t),t)=D||x,x=1..NAE)}:

Dvars2:={seq(rhs(Dvars1[x])=lhs(Dvars1[x]),x=1..NAE)}:

icsn:=seq(subs(y||x(0)=y||x(t),icodes[x]),x=1..NODE),seq(subs(z||x(0)=z||x(t),icaes[x]),x=1..NAE):

for j from 1 to NAE do

EQAEX||j:=subs(Dvars1,eqode,icsn,Dvars2,lhs(EQAE||j))=rhs(EQAE||j):

end do:

Sys:={seq(EQODE||x,x=1..NODE),seq(EQAEX||x,x=1..NAE),seq(icodes[x],x=1..NODE),seq(icaes[x],x=1..NAE)}:

if Xexplicit=1 then

sol:=dsolve(Sys,numeric,maxfun=0):

else

sol:=dsolve(Sys,numeric,stiff=true,implicit=true,initstep=0.1,compile=true,maxfun=0):

end if:

 

for XX from 1 to NODE do

a||XX:=odeplot(sol,[t,y||(XX+1)(t)],1..subs(pars,tf),color=red): end do:

for XX from NODE+1 to NODE+NAE do

a||XX:=odeplot(sol,[t,z||(XX-NODE)(t)],1..subs(pars,tf),color=blue): end do:

display(seq(a||x,x=1..NODE),a||(NODE+NAE-1),a||(NODE+NAE),axes=boxed);

 

E

 

Download example4improved.mws

@acer 

 

works like a charm in Maple2015. Rosenbrock methods keep the same coefficient matrix at different stages, and this is extremely useful for that.

@roman_pearce 

 

Roman, I will get back to you on this soon after I gather references. Meanwhile Maple 14 is able to do sparse Linear solve for Digits: = 20 (software float), but Maple 2015 cannot do Linear Solve for more than 15 Digits using the sparse option. Is this because Maple moved to a different sparse solver? It will be great to continue providing high precision solvers.

Thanks

 

@ 

According to this post from roman_pearce it can't be done. Perhaps Maple will consider providing this option in the future releases

@acer 

Stiff ODE solvers (Rosenbrock) involve using the same coefficient matrix with different rhs. For example AX=b, followed by AX1=b1 (where b1 = linear function of b), etc. Simple LU deomposition will allow for one LU decomposition followed by multiple backsubstiution. Is it possible to do this with SparseLU or SparseDirect option?

Thanks

 

Being a professor, it pains me to see non-deterministic algorithms (Genetic algorithms, etc) being pushed everywhere. I have seen students directly trying stochastic approaches for problems that are well defined and are easily solved using gradient based techniques.

There is a place for non-gradient approaches, but just using stochastic, genetic algorithms and calling them global optimzation is wrong.

This is unfortunately true even with GlobalOptimization package provided by Maple as well.

 

@Greg3141 

 

that will show the content of the code and the steps involved (in regular maple).  I don't think one can hide the content in Maple. This is a serious problem when confidential/classified information are to be hidden while sharing the codes.

 

@Carl Love 

 

This is how Maple solves DAEs. (Differentiate algebraic equations, arrive at dy/dt = f(y)) and then use AE as a residue. 

@Carl Love 

 

That has been my biggest complaint. I am ok if Maple is not able to sovle implicit DAEs, but it says it solves DAEs. But it only solves DAEs by converting DAEs to ODEs (unless the Algebraic equation is simple linear equation or linearly solvable in terms of ODE variables).

 

Write any DAE and solve with any DAE solver in Maple. Use infolevel = 10, you will see YP for every variable including algebraic.

 

That is the sad part (incorrect help file)

 

 

@ patient

YP[2] = (-6*X*Y[3]*Y[1]^2*Y[4]^2+Y[1]*Y[2]*(-3*X^2*Y[1]^(1/2)+Y[4]^2*Y[3]*(X^2-8*Y[4]^2))+6*Y[4]^2*X*Y[2]^2*(Y[4]^2*Y[3]-Y[1]^(1/2)))/(Y[1]^(1/2)*X*(4*Y[4]^4*Y[1]^(1/2)*Y[3]+X^2*Y[1]+2*Y[4]^2*X*Y[2]))

YP[3] = Y[3]*(-3*X*Y[1]^2+2*Y[4]^2*(-Y[4]^2*Y[1]^(1/2)*Y[3]+X*Y[2]+Y[1])*Y[2])/(Y[1]*(4*Y[4]^4*Y[1]^(1/2)*Y[3]+X^2*Y[1]+2*Y[4]^2*X*Y[2]))

YP[1] = Y[2]

 

The variables are S, S' and U. Y[4] is x0.

X is the independent variable.

@tomleslie 

Maple does not have a DAE sovler. If Maple solves with rkf45dae it probably means that is solves dy/dt =f with algebraic equations may be as a constraint.

 

Patient, convert all of your equations to dy/dt =f . If needed differentiate the algebraic equations, then you will be close to what Maple does.

 

If some one posts *.mws worksheet that works, I will try. ( or at least in that format, 1D math)

 

@momoklala 

 

if code works for B = -0.5, and L = particular value, use that as approxsoln for B =0. 

 

 

@Alejandro Jakubi 

 

Acer and Alejandro, I kind of agree partially with both of you regarding openness of Maple. I started using Maple from Maple IV onwards as a student, and have slowly started learning about analytical methods/numerical methods, odes and DAEs. The kind of problems I face in my research even now cannot be solved by most solvers/platforms all the time. But Maple has helped me learn different nuances in numerical methods easily ( I couldn't handle mathematica which was not userfriendly for me).

While Maple is reasonably open, (Explore command is great from LibraryTools), documentation is bad and many things are not properly explained. To me, it looks like the focus right now is to get the attention of kids doing middle school or high school math, with no attention/care towards users having to solve large number of equations/ODEs/DAEs.

Another thing I hate is the inaccuracy in the claims. Sorry to repeat again, while Maple claims it can solve DAEs, it can only solve ODEs and should acknowledge the same. The worst thing is users reading about Jeff Cash who developed mebdfi code might wrongly assume that he is wrong and his method cannot solve DAEs. This is not true, it is just the Maple has implemented it wrongly. Ironically, Maplesim can solve DAEs. Even if you have license for both, you can't use Maplesim's DAE solver without creating a custom component which will take an hour or  more to create. For a change instead of just complaining, you will see me posting solutions for DAE problems in Maple in the next few days or weeks (just now a research paper was accepted, this approach is dfifferent from Maplesim or any other DAE solvers published as of today).I don't know where else Maple has made such mistakes/inaccurate claims.

Recently Acer and myself arrived at a compiled Newton Raphson. I can't understand why Maple cannot include such a thing in its fsolve (which consumes too much RAM).

FWIW, my view on where Maple can play a big role and is not doing it.

Instead of Codegen being able to convert Maple code to C, it should do it the otherway around. (Take an external C file and convert to Maple. Most algorithms have specified number of input variables, parameters and expected output variables). I know that there is a way to do externalcalling, but it is not well documented and difficult. Why do I want this? This could mean that Maple does not need to write anything new, just use existing good or new things coming out in C or Fortran directly.

Maplesim web based approach is a good step. You want to be able to share your results with others who may not have access to the software. A good business model should allow this. Only then potential users will sign up for license. If having Maple is a prerequisite to share results, then how can one expand the user pool?

 

 

First 11 12 13 14 15 16 17 Page 13 of 17