Venkat Subramanian

296 Reputation

13 Badges

11 years, 329 days

MaplePrimes Activity

These are replies submitted by

z goes only from - 1 to 0. 

u1,theta1,N solved from -1 to 0.
u2, theta2 are also solved from -1 to 0 (for z)
I think you have 3 variables in region 1 and 2 in region 2. variables and fluxes are continuous at y = 0. 

it will be great in the future if you state all the bcs at y = -1 first, then at y =0, and then at y = 1, for the convenience of coding. 

Until Maple 17 or so D(u[1])(0) worked. In the later version, you have to use u1.


BVP is solved in z. z varies from -1 to 0.

eq1:=diff(u1(z),z,z) + 1/2*diff(N(z),z) + 1/2*theta1(z) = 0;

eq1 := diff(u1(z), `$`(z, 2))+(1/2)*(diff(N(z), z))+(1/2)*theta1(z) = 0


eq2:=diff(N(z),z$2) - 2/3 * (2*N(z) + diff(u1(z),z)) = 0;

eq2 := diff(N(z), `$`(z, 2))-(4/3)*N(z)-(2/3)*(diff(u1(z), z)) = 0


eq3:=diff(theta1(z),z$2) = 0;

eq3 := diff(theta1(z), `$`(z, 2)) = 0


eq4:=diff(u2(z),z$2) + theta2(z) = 0;

eq4 := diff(u2(z), `$`(z, 2))+theta2(z) = 0


eq5:=diff(theta2(z),z$2) = 0;

eq5 := diff(theta2(z), `$`(z, 2)) = 0



bcs := u1(-1) = 0, theta1(-1) = 1, N(-1) = 0, u2(0) = 0, theta2(0) = 0, u1(0) = u2(-1), (D(N))(0) = 0, theta1(0) = theta2(-1), (D(theta1))(0) = (D(theta2))(-1), (D(u1))(0)+(1/2)*N(0) = (1/2)*(D(u2))(-1)




[z = -1., N(z) = 0., diff(N(z), z) = -0.904317007826591e-1, theta1(z) = 1.00000000000000, diff(theta1(z), z) = -.500000000000000, theta2(z) = .500000000000000, diff(theta2(z), z) = -.500000000000000, u1(z) = 0., diff(u1(z), z) = .371898088450732, u2(z) = .172870489765202, diff(u2(z), z) = -0.620382309853562e-2]



[z = 0., N(z) = -0.228483654680803e-1, diff(N(z), z) = 0., theta1(z) = .500000000000000, diff(theta1(z), z) = -.500000000000000, theta2(z) = 0., diff(theta2(z), z) = -.500000000000000, u1(z) = .172870489765202, diff(u1(z), z) = 0.832227118477233e-2, u2(z) = 0., diff(u2(z), z) = -.256203823098536]







Please confirm that this is ok.

Download threepointBVPmodify.mws


Since this example has only two inputs, f converts the vector to scalars.

It will be great to have fdiff which takes obj which takes a vector and performs the calculation. When we have 10 or 20 inputs to the objective, it is a pain to convert the model to scalars from vectors. In that case, why can't we write obj in the operator form and use it directly?



Please note that stiff error comes when stiff=true is used. This is the failure of Rosenbrock methods, probably from Maple's failure to find the analytic Jacobian.

For the OP, you may be better off writing your system as a system of first order equations (Maple does this before going to dsolve). Also try implicit =true to see if that helps. This can be done only for first order ODEs.


Thanks for your solution. I believe the NLPSolve with method=nonlinearsimplex has no issues (it can't deal with bounds and constraints). It would be great if NLPSolve can use an inbuilt wrapper to make the optimization work for simple returns from the procedures without these _passed numeric returns. I believe the issue may be due to the calls for fdiff. 
Please note that GlobalOptimzation package calls the sqp optimizer, and has no issues with procedures without the extra _passed commands. 

@ I am able to run classicworksheet with maple2021

I had to go to to folder

run cwmaple.exe

This is actually better than previous versions where I could only run 32bit version in classic worksheet mode. The help file says that cwmaple.exe exists inside bin.x86 folder,  but it is in folder


I have classic worksheet working for Maple2019 (a little unstable, but it works). So classic worksheet (32bit) was dropped either in 2020 or 2021.

Unfortunately no more classic worksheet. I understand that there may not be a demand for 32 bit installations, but I am a big fan of classic worksheet interface, please add it if possible.


Please note that this statement (IVP numerical methods are more stable compared to BVP methods) is incorrect. Actually, BVP methods are more stable to solve a given BVP problem. It is easy to show many examples to prove the same. For example, 

take dy1/dt = y2

dy2/dt = tau^2*y1 + (Pi^2-tau^2)*sin(Pi*t)

IVP: t=0, y1= 0, y2 = Pi

BVP: t= 0, y1 = 0; t = 1, y1 = 0

(Credit: Dr. Biegler at CMU)

This model is very stable with a BVP approach, but not with IVP approach. Shooting method is not useful for largescale or stiff problems. At least multiple shooting is needed.

Add some noise to ics, constants, you will see BVP methods behaving much better.



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



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


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




Digits := 15



N := 5

h := .166666666666667



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:








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.


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.




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.


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