Solving DAEs in Maple

As I had mentioned in many posts, Maple cannot solve nonlinear DAEs. As of today (Maple 2015), given a system of index 1 DAE dy/dt = f(y,z); 0 = g(y,z), Maple “extends” g(y,z) to get dz/dt = g1(y,z). So, any given index 1 DAE is converted to a system of ODEs dy/dt = f, dz/dt = g1 with the constraint g = 0, before it solves. This is true for all the solvers in Maple despite the wrong claims in the help files. It is also true for MEBDFI, the FORTRAN implementation of which actually solves index 2 DAEs directly. In addition, the initial condition for the algebraic variable has to be consistent. The problem with using fsolve is that one cannot specify tolerance. Often times one has to solve DAEs at lower tolerances (1e-3 or 1e-4), not 1e-15. In addition, one cannot use compile =true for high efficiency.

The approach in Maple fails for many DAEs of practical interest. In addition, this can give wrong answers. For example, dy/dt = z, y^2+z^2-1=0, with y(0)=0 and z(0)=1 has a meaningful solution only from 0 to Pi/2, Maple’s dsolve/numeric will convert this to dy/dt = z and dz/dt = -y and integrate in time for ever. This is wrong as at t = Pi/2, the system becomes index 2 DAE and there is more than 1 acceptable solution.

We just recently got our paper accepted that helps Maple's dsolve and many ODE solvers in other languages handle DAEs directly. The approach is rather simple, the index 1 DAE is perturbed as dy/dt = f.S ; -mu.diff(g,t) = g. The switch function is a tanh function which is zero till t = tinit (initialization time). Mu is chosen to be a small number. In addition, lhs of the perturbed equation is simplified using approximate initial guesses as Maple cannot handle non-constant Mass matrix. The paper linked below gives below more details.

http://depts.washington.edu/maple/pubs/U_Apprach_FULL_DRAFT.pdf  

Next, I discuss different examples (code attached).

Example 1: Simple DAE (one ODE and one AE), the proposed approach helps solve this system efficiently without knowing the exact initial condition.

Hint: The code is given with a semicolon at the end of the most of the statements for educational purposes for new users. Using semicolon after dsolve numeric in classic worksheet crashes the code (Maplesoft folks couldn’t fix this despite my request).

Example 2:

This is a nickel battery model (one ODE and one AE). This fails with many solvers unless exact IC is given. The proposed approach works well. In particular, stiff=true, implicit=true option is found to be very efficient. The code given in example 1 can be used to solve example 2 or any other DAEs by just entering ODEs, AEs, ICs in proper places.

Example 3:

This is a nonlinear implicit ODE (posted in Mapleprimes earlier by joha, (http://www.mapleprimes.com/questions/203096-Solving-Nonlinear-ODE#answer211682 ). This can be converted to index 1 DAE and solved using the proposed approach.

Example 4:

This example was posted earlier by me in Mapleprimes (http://www.mapleprimes.com/posts/149877-ODEs-PDEs-And-Special-Functions) . Don’t try to solve this in Maple using its dsolve numeric solver for N greater than 5 directly. The proposed approach can handle this well. Tuning the perturbation parameters and using compile =true will help in solving this system more efficiently.

Finally example 1 is solved for different perturbation parameters to show how one can arrive at convergence. Perhaps more sophisticated users and Maplesoft folks can enable this as automatically tuned parameters in dsolve/numeric.

Note:

This post should not be viewed as just being critical on Maple. I have wasted a lot of time assuming that Maple does what it claims in its help file. This post is to bring awareness about the difficulty in dealing with DAEs. It is perfectly fine to not have a DAE solver, but when literature or standard methods are cited/claimed, they should be implemented in proper form. I will forever remain a loyal Maple user as it has enabled me to learn different topics efficiently. Note that Maplesim can solve DAEs, but it is a pain to create a Maplesim model/project just for solving a DAE. Also, events is a pain with Maplesim. The reference to Mapleprimes links are missing in the article, but will be added before the final printing stage. The ability of Maple to find analytical Jacobian helps in developing many robust ODE and DAE solvers and I hope to post my own approaches that will solve more complicated systems.

I strongly encourage testing of the proposed approach and implementation for various educational/research purposes by various users. The proposed approach enables solving lithium-ion and other battery/electrochemical storage models accurately in a robust manner. A disclosure has been filed with the University of Washington to apply for a provisional patent for battery models and Battery Management System for transportation, storage and other applications because of the current commercial interest in this topic (for batteries). In particular, use of this single step avoids intialization issues/(no need to initialize separately) for parameter estimation, state estimation or optimal control of battery models.

Download SolvingDAEsinMaple.mws

This is a very simple suggestion that weights heavily on my enjoyment of Maple.

I'm not sure if it is only me and my students but it is really tricky to resize a code edit region. Trying to get a hold of the contour in an exercise in patience! Can anyone fix this?

Thanks!

Stéphane

I do not think the current API for dsolve when asking for series solution is done right. If one wants to obtain a series solution for an ODE, but wants different order than the default 6, now one must set this value using a global setting before making the call, like this:

eq:=diff(y(x),x$2)+y(x)=0;
Order:=10;
dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series');

It would be better if options to calls are passed along with the other parameters in the call itself. Something like

dsolve({eq,y(0)=1,D(y)(0)=0},y(x),type='series',order=10);

This should also apply to any command that takes Order, such as series(cos(x),x=0,order=10);

Passing options and values to functions using global and environment variables is not safe and not a good way to go about it as it can cause programming errors.

@ecterrab 

I figured I'd start a new thread for odd things I come across whilst using the new physics package. 

I have found this, and am not sure if it is expected. 

Download Simplify2.mw

Did Maplesoft make a statistics about questions on MaplePrimes?

There must be some commands frequently meet problems.

I listed 10 commands or packages need to be enhanced much more.

1.fsolve.

eg, roots of transcendental equation on complex plane & system of nonlinear equations.

2.int.

Both numerical integration & symbolic integration need further developed.

3.Eigenvalues&LinearSolve.

large scale computing, everyone needs it.

4.dsolve&pdsolve.

dsolve needs to be enhanced for solve more equations in parallel.

As mathematica 10 contains the FEM package, I think Maple needs to do something, develop a new package needs more time than make the existing pdsolve command better.

5.Threads package.

parallel programing is the future. BTW, why don't make the task model much easier to use, like the mathematica, just a command?

6.plot&plot3d

there is much space to go further about the two commands.

7.Optimization package or Statistics package.

In fact, both packages need to be further enhanced, especially for Optimization package.

Although some modern algorithms about optimization&statistics is much easier to see on the internet, the

Maplesoft should not stop, but go further, add them in Maple.

Above is just my opinion.

Thanks for your attention.

The problem is when initializing a Matrix with a list of strings. The worksheet excerpt below shows the normal behavior using a list of integers to initialize a square matrix: the successive list elements fill the matrix by rows.

Then trying the exact same thing with a list of strings instead of integers gives an error message!
This is not right. While it is an odd and likely rare problem, it would be better fixed.

x := [i $ i=1..25];
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25]
Matrix(5,5,x);
 

y := map(convert,x,string);
["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25"]
Matrix(5,5,y);
Error, (in Matrix) initializer defines more columns (25) than column dimension parameter specifies (5)

Greetings to all.

I would like to share a brief observation concerning my experiences with the Euler-Maclaurin summation routine in Maple 17 (X86 64 LINUX). The following Math StackExchange Link shows how to compute a certain Euler-MacLaurin type asymptotic expansion using highly unorthodox divergent series summation techniques. The result that was obtained matches the output from eulermac which is definitely good to know. What follows is the output from said routine.

> eulermac(1/(1+k/n),k=0..n,18);
     1       929569        3202291        691                O(1)
O(- ---) - ----------- + ----------- - --------- + 1/1048576 ----
     19             15            17          11              19
    n      2097152 n     1048576 n     32768 n               n

                                           n
                                          /
        174611      5461        31       |      1           17        1
     - -------- + --------- + ------- +  |   ------- dk - ------- + ------
             19          13         9    |   1 + k/n            7        5
       6600 n     65536 n     4096 n    /                 4096 n    256 n
                                          0

         1       1
     - ------ + ---- + 3/4
            3   16 n
       128 n

While I realize that this is good enough for most purposes I have two minor issues.

• One could certainly evaluate the integral without leaving it to the user to force evaluation with the AllSolutions option. One can and should make use of what is known about n and k. In particular one can check whether there are singularities on the integration path because we know the range of k/n.
• Why are there two order terms for the order of the remainder term? There should be at most one and a coefficient times an O(1) term makes little sense as the coefficient would be absorbed.

You might want to fix these so that the output looks a bit more professional which does enter into play when potential future users decide on what CAS to commit to. Other than that it is a very useful routine even for certain harmonic sum computations where one can use Euler-Maclaurin to verify results.

Best regards,

Marko Riedel

Maple 18, document mode. Please try this:

>restart;
>some_long_name:=3:
>f(som)

Now, as the cursor is right after the "m" above, hiting the ESC key to autocomplete, nothing happens. Only a beep, meaning Maple does not know about the variable name I wanted to expand to, which is "some_long_name". What I do now is manually add a space by  pushing the ")" away, like this:

  f(som )

And now put the mouse back in front of the "m" and now hit ESC. Maple now can see the name I wanted.

Really? This is very bad design. You might ask, why did I close the () first? Well, I like to start by writing () then fill it in, or it can be that I am changing things, and wanted to write a new name, and the () was allready there. I use name autocomplete alot, and it is very annoying to not be able to use it if happens that there is a ")" or "]" right next to the name.

Having to keep pushing the ")" away, so that Maple can see the name makes no sense since ")" can't part of a name.

In Mathematica it just works. I can do

someLongName = 3;

Then in a new cell,

f[som] and hit CTRL-K  when mouse at the letter "m" just typed it in, and it will autocomplete.  I did not have to push "]" away first.

Please Maple, fix your user interface so it is less awkward and annoying to use. This is version 18 of the software, and not version 1.

A little late in the game, with no Maple 18 wishlist yet I'd like to get it started

granted all wishes for previous versions not yet applied should still, hopefully, be under consideration.  

1 - The abitlity to angle x-axis labels.

2 - Textplot option for text to rotate with graph as it is rotated.  Currently all textplot texts are held to the horizontal

It would be nice if Maple had a procedure which could turn a procedurelist into a listprocedure.
I use these words in the sense they are used in dsolve/numeric.
Thus by a procedurelist I mean a procedure returning lists (of numbers).
By a listprocedure I mean a list of procedures all having the same formal parameters and each returning one number.
Thus a `convert/listprocedure`should accept a procedurelist p as input and give as output the corresponding listprocedure [p1,p2, ... pN] where N is the number of elements in the output from p.

I tried making a `convert/listprocedure` myself and in doing so found that it was not totally trivial.
I had lots of problems but did end up with something that seems to work.

convert-listprocedur.mw

But my main point is that Maple ought to have some such facility either as described available to the user or by changing fsolve, complexplot or what have you, so that procedurelists are accepted.

I really very much like this package since it works just perfectly converting units and also supports CGS for people doing calculations in theoretical physics. I have just a suggestion concerning the formattig: As is explained in many physics books and also at the units standard website http://physics.nist.gov/, units shouldn't be placed in square brackets since these brackets are defined such that [m] = kg means: unit of m is kg. To write is m = 3 [kg] is a very common and unharmful mistake, also used very often when labelling axes in figures, e.g. m [kg] instead of the correct m/kg or m (kg). Another issue is that units shouldn't be printed italic but upright and any physics journal requests its authors to consider this rule. Again, this is not a "mistake" but on the other hand the Maple pretty-print output is already so nice that considering such rules would even make it perfect.

I am fully aware that omitting the brackets could be considered harmful if someone has an algebraic expression in front of the unit where symbols and units could be mixed. On the other hand having the usual interspace between unit and expression, the unit being upright and maybe even coloured (analogously to e.g. operators in the physics package) units and symbols could be easily held apart.

Bug in type(HFloat(-infinity),  pos_infinity). Negative infinity is incorrectly recognized as positive one:

s:=HFloat(-infinity);
                   HFloat(-infinity)

type(s, neg_infinity);
type(s, pos_infinity); # bug
                    ...

Greetings to all.

I have been using the numtheory package for quite some time now and it has helped me advance on a number of problems. Recently an issue came to my attention that I have known about for a long time but somehow never realized that it can be fixed. This is the fact that the numtheory package does not know about Dirichlet series, finite and infinite. Here are two links:

Greetings to all. I will keep this brief and to the point. I would like to point out a certain deficit in the integral transform package. I have recently been calculating some Mellin transforms at this link and the base functions are of the following type.

g := (p, q) -> 1/(x+p)/(x+(p*q-1)/q);

Now to see the deficit here are some Mellin transforms that...

I am pretty new to programming in Maple and trying to learn the language. One thing I needed was how to handle lists in a specific way. I have used many hours on this issue now and even asked questions in the Questions section. Other users helped me with clever tricks and I am grateful for that. What is left is however a feeling that one need to know a lot of rather Ad Hoc commands to accomplish your tasks. I looked for operations in the ListTools package, going through all...