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

@Preben Alsholm 

I am not surprised. I am actually pleased that I had the patience to browse through a 2D math input code and manage a working code out of it. Thanks for the catch.

I will work on this code more later. I think there is a singularity at both eta = 0 and eta = N. midrich method avoids the end points. Exp(1/0) should not occur in iteration, but the nearest node point probably has issues.

For problems like this providing a mesh (0,1e-6...0.2, 0.4, 0.99, 0.999999, 1)total length might work.

@acer 

used evalf for creating the F procedure (Sometimes users might provide 2x-x^2 instead of 2.x-x^2)

Also, as you are not adding the error, you don't need to divide by N in the Newton procedure. This should be faster

Let us get greedy next,

For a given list of equations and variables with guess , eqs, ICs, can we create a dll or downloadable package/library addition for this?

Download newtoncompiled6.mw

(1) First please use only Maple input and 1D math if you can. Classic worksheet is the easiest and least buggy. I spent a lot of time on this problem and most of which was just spent on copying and pasting your equations. (Is there a way to take a mw worksheet or 2D math code to automatically convert to 1D math input and classic worksheet?). I wish common sense prevailed and Maple never moved to mw.

I attach a partial solution. I think this is a semi-infinite domain problem. I am using N = 2, 10 may be an overkill.

Also, at N =2, physical intuition suggests that flux of theta can be taken to zero (I am assuming that you are solving for velocity and temp).

Low values of epsilon (that is a crazy symbol, I never used it and it was a pain to copy and paste), gives a solution which can be used as an approxsoln(guess) for large values. I leave the rest to you.

If you want solution for your original bc, take my solution and use as approxsoln for your case. might work. Another trick is to write continuation on theta(N) based on my bc.

Few tips, increase Digits, use continuation parameter based on the physics. (Ideally a linear problem or a simple problem from which you can grow to your problem)

Download VT1.mw

@ 

At least in Maple 14, you will see that fsolve consumes huge amount of RAM compared to the discussed approach. I think this is because fsolve stores the zeros of the Jacobian.

@acer 

If you manually copy paste F , Jac and s3 and Newton in a single procedure and compile, that will be faster. (I did this before, can't find it).

I think it is interesting to note that both Newton (original code) and Newtonc have similar speed for a single iteration, but I expect that a single compiled procedure (without global calls) will be faster.

@acer 

Yes, I typically remove 0 for compiling otherwise it is too slow (and RAM intensive for large problems). While this is good, unfortunately, this approach is not useful for dsolve/numeric/compile for DAEs.  Please check my PM for DAE, as  I want this thread focussed strictly on Newton Raphson.

Also, fsolve does not provide option for tolerance. So, providing tolerance as 1e-4 or 1e-6 is more than enough for many problems. That is much appreciated.

You don't need to do s3 or s3c inside a dsolve/numeric. You can easialy call 'dsolve/numeric/SC/Decomp' followed by 'dsolve/numeric/Solve', both of which are recognized in dsolve/numeric.