Venkat Subramanian

276 Reputation

13 Badges

10 years, 212 days

MaplePrimes Activity


These are replies submitted by

@janhardo 

I meant a tablet/laptop not the browser.

If you install 32 bit Maple, you should be able to use classic worksheet mode. The stability is anyone's guess. But worth trying

@nm 

I think when Maple was started, the goal was perhaps to keep the source codes transparent. New additions (maplesim, stochastic packages, etc) have hidden codes that are nearly impossible to dig. 

Mathematica has everything hidden. MATLAB has many codes that are transparent, but not FFT. Perhaps true for A\b as well.

@acer 

 

I tried these commands, but I tried in Maple14. I see that this works now in later versions of Maple. thanks.

 


I don't know the physics of your model. See a simple ODE system with z(t) as a possible algebraic variable.

If z(t) =cos(Pi/2t), the code will hit singularity as the coefficients of time derivatives might go to zero.

First z(t) is explicitly specified and then I solve the same model with z(t) as an algebraic variable. Doing this type of analysis will help you understand if any of the coefficients are going to zero or nearly zero causing possible singularity in your codes.

You might want to try stricter tolerance in the code Preben sent to see if that helps as well. For this case z(0)=2, but Maple is able to fix the mistake even if we give a wrong guess for z(0). For nonlinear problems, you have to provide correct initial condition for z.

 

restart;

z(t):=piecewise(t<5,1+cos(Pi*t/2),1);

"z(t):={[[1+cos((Pi t)/(2)),t<5],[1,otherwise]]"

(1)

#z(t):=cos(Pi*t/2);Use of z = cos(Pi*t/2) will make it singular after sometime

eq1:=(z(t))*diff(y1(t),t)+(1+z(t))*diff(y2(t),t)=exp(-t)*sin(t);

eq1 := piecewise(t < 5, 1+cos((1/2)*Pi*t), 1)*(diff(y1(t), t))+(1+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1))*(diff(y2(t), t)) = exp(-t)*sin(t)

(2)

eq2:=(2+z(t))*diff(y1(t),t)+z(t)*diff(y2(t),t)=t*sin(t);

eq2 := (2+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1))*(diff(y1(t), t))+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1)*(diff(y2(t), t)) = t*sin(t)

(3)

 

sol:=dsolve({eq1,eq2,y1(0)=1,y2(0)=0.5},type=numeric,implicit=true,maxfun=0,stiff=true):

sol(0);sol(1);

[t = 0., y1(t) = 1., y2(t) = .500000000000000]

[t = 1., y1(t) = 1.06023325551433, y2(t) = .565134897638036]

(4)

plots:-odeplot(sol,[t,y1(t)],0..15,thickness=3,axes=boxed);

 

plots:-odeplot(sol,[t,y2(t)],0..15,thickness=3,axes=boxed);

 

 Next solve as a DAE

z(t):='z(t)';

z(t) := z(t)

(5)

eq1;eq2;

piecewise(t < 5, 1+cos((1/2)*Pi*t), 1)*(diff(y1(t), t))+(1+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1))*(diff(y2(t), t)) = exp(-t)*sin(t)

(2+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1))*(diff(y1(t), t))+piecewise(t < 5, 1+cos((1/2)*Pi*t), 1)*(diff(y2(t), t)) = t*sin(t)

(6)

restart;

eq3:=z(t)=piecewise(t<5,1+cos(Pi*t/2),1);

eq3 := z(t) = piecewise(t < 5, 1+cos((1/2)*Pi*t), 1)

(7)

eq1:=(z(t))*diff(y1(t),t)+(1+z(t))*diff(y2(t),t)=exp(-t)*sin(t);

eq1 := z(t)*(diff(y1(t), t))+(1+z(t))*(diff(y2(t), t)) = exp(-t)*sin(t)

(8)

eq2:=(2+z(t))*diff(y1(t),t)+z(t)*diff(y2(t),t)=t*sin(t);

eq2 := (2+z(t))*(diff(y1(t), t))+z(t)*(diff(y2(t), t)) = t*sin(t)

(9)

sol:=dsolve({eq1,eq2,eq3,y1(0)=1,y2(0)=0.5,z(0)=1},type=numeric,maxfun=0,stiff=true):

sol(0);sol(1);

[t = 0., y1(t) = 1., y2(t) = .500000000000000, z(t) = 2.]

[t = 1., y1(t) = 1.06023333732083, y2(t) = .565134673122847, z(t) = .999999999931483]

(10)

plots:-odeplot(sol,[t,y1(t)],0..15,thickness=3);

 

plots:-odeplot(sol,[t,y2(t)],0..15,thickness=3);

 

plots:-odeplot(sol,[t,z(t)],0..15,thickness=3);

 

 


 

Download DAEapproach.mws

Your model has only 6 ODEs, but complicated rhs. Try to define Lm, id and everything as additional algebraic variables.

Then try stiff=true, implicit=true,compile=true (try stiff=true, compile=true, optimize =fast first to see if that helps).

If you post the code in mws form (classic worksheet) mode without using any procedures on the rhs, I can try to do this and check.

There is a catch here as well as Maple doesn't have a direct DAE solver so it may not work in particular for piecewise functions.

Don't try to solve ODEs to get dy/dt = f, that just makes it worse in this case. 

@acer 

I agree with you on the pedagogy aspects. Is there any gain in CPU time, memory or anything else by going to embedded coding?

@ looks like fsolve was updated in Maple 2018. I am yet to test it. 

I have used most of them and I can provide some thoughts.

@eithne 

At this point in time, it seems to me that 2019.2 should be withdrawn. There are some benefits, but negatives are significant.

 

@acer 

Most users would know a good initial guess and Globalsolve should attempt to call a localsearch based on the initial guess provided by the user first. This is not just with Maple, many packages seem to ignore the initial guess provided by the user for the GA.

@acer 

Any global solver should check results from a local solver before displaying the final answer. The word Global is abused and there are no guarantees for global optima, in particular for blackbox objectives.

This may be because the global optima solver probably ignores the initial guess provided and starts searching. It should ideally first do a gradient-based search based on the initial guess provided. 

But we can't win against the craze/push for GA.

Though mapleprimes will help you with maple related questions, I will be surprised if you find help for this. In my opinion, Maple's initial goal was to call define_external and build on external routines. They should have just done this instead of rewriting (sometimes wrong) codes for dsolve/numeric, pde/numeric, etc. Perhaps, compiler changes make this complicated, but many commands won't work well for define external (for fortran or C). Though I like maplesim, in my opinion, Maple folks dropped the ball on this after Maplesim came online.

For example, I would like to see RK, or Fortran codes (MEBDF, eg.,) in Fortran and someone explaining how this gets converted to dll files in Maple. If we get help for this (including sparse matrices/vectors), Maple can get very powerful (by calling all possible open-source ODE/DAE/PDE/FEM solvers, optimizers).

Edits: I have learned a lot from the current implementation of ODE/DAE solvers in Maple and the transparent nature of implementation helps in identifying errors when the solvers fail. The same thing can't be said about NLPSolve. It frequently fails for gradient-based inequality constrained optimizations for blackbox objectives. If the plan is to call external routines, implementing them correctly is very important.

@mmcdara 

I had been wanting to write to this post.

We should write in flux form and assign conductivity to the harmonic mean of values at two neighbouring nodes. This applies in 2D as well.

 

@tomleslie 

I believe a procedure is the best way to approach this. But there are situations when the OP's approach is better. For example, let us assume that we need to call fsolve or Roots or something similar and it fails. In the OP's approach, it will be easier to debug, and in a procedure, it will be a pain.

Once someone is absolutely sure of all the commands in a code, then the procedure makes sense. The OP may not have the ability to share the code (copyrights/etc).

@Daniel Skoog 

Please don't stop classic worksheet. Though it is unstable now, that is what I use. For me it is 100x faster than any of the Java based interfaces.

Also, it will be great to have a fsolve that works and gives answers for a specified tolerance. Provide an option that uses the Newton Raphson for mutliple variables which is guaranteed to give a better answer compared to the starting initial guess just with Digits:=15 without increasing Digits to seek unrealistic tolerances. 

 

1 2 3 4 5 6 7 Last Page 1 of 14