Venkat Subramanian

386 Reputation

13 Badges

15 years, 243 days

MaplePrimes Activity


These are replies submitted by

@acer 

The current linear sparse direct solver is good. If we can call that in evalhf or compiled procedure (making sure of arrays, etc), that will be of great help for solving large scale problems. Creating an interface to MUMPS or PARDISO will be a great step up (to enable the use of all the cores in the computer).

 

@acer 

I haven't tested mathematica for parallel Sparse Solvers. But MATLAB has this.

@acer 

By arbitrarily defining a 2D grid for linear elliptic problems, I don't see SparseDirect solvers using more than one core in my computer with any version of Maple. Can you please confirm if you are able to run Maple for the following simple script (for 1D problem) for very large N using all the cores? Try N =999, 999999, etc


 

restart;

Digits:=15;

Digits := 15

(1)

N:=5;h:=1.0/(N+1);

N := 5

h := .166666666666667

(2)

A:=Matrix(1..N,1..N,datatype=float[8],storage=sparse):

for i from 1 to N do A[i,i]:=-2/h^2-1.0:
if i >1 then A[i,i-1]:=1/h^2:end:
if i<N then A[i,i+1]:=1/h^2:end:od:

b:=Vector(1..N,storage=sparse,datatype=float[8]):

b[1]:=-1/h^2:

y:=LinearAlgebra:-LinearSolve(A,b,method=SparseDirect):

y[N/2+1/2];

.443527611152833

(3)


 

Download LargeN.mws

 

1. Maple doesn't have a DAE solver, the help file should state the same and instead state that DAEs are converted to explicit ODEs and solved. If possible add the DAE solver in Maplesim to Maple.

2. Make NLPSolve:-Optimize more robust and accurate. For example, if you call a procedure with 2 arguments or more and optimize with bounds and constraints, it will inevitably say -first-order conditions met and optimization not continued. This seems to be a bug in either numerical Jacobian calculation or unfortunately a trick to force users to buy the Global Optimization package. I hope that it is not the case.

3. Update ODE, DAE, BVP, and PDE numerical solvers with sparse linear algebra.

4. Enable parallel computing for sparse linear solver (call MUMPS or PARDISO).

5. Make 1D input the default or at least provide that option when installing the software or opening Maple the first time.

@ecterrab 

Liked your post. The goal should be addressing a problem. 

I would like to see 1D as the default mode as opposed to 2D for Maple? Or at least ask the user during installation to choose 1D as opposed to 2D.

In my humble opinion, Maple should not have done anything beyond the classic Maple worksheet in terms of interface. When that consultant went out of business, Maple should have invested in the same. Typing in Maple 8 or earlier versions was just like typing in a text file. It consumed very little RAM and was very responsive.

 

 

@Aminaple 

take my example 1 and modify it to your example or problem. If needed define your model as first order equations (convert second order to first order equations)

The solver uses implicit midpoint approach for first-order BVPs. (Second-order or higher-order ODEs are converted to first-order systems). The algorithm is given below.

dy/dx = f

y1= y0+h*f(ymid) where ymid is (y1+y0)/2, h is the element size. h = 1/N, N the number of elements. 

Caveat: The Richardson extrapolation method assumes second order of accuracy, but if your model is stiff or ill-conditioned or nearly singular, this will sometimes give wrong/sub-optimal results as the order of accuracy is 1 and the solver does not take care of this.

 

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

2 3 4 5 6 7 8 Last Page 4 of 17