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

I'd like to pay attention to Maple packages for symmetric functions, posets, root systems, and finite Coxeter groups by John Stembridge. That was estimated as well qualified in Mathoverflow.

Dr Lopez,

Thank you for your very interesting webinar on eigenpairs today. I especially valued your "aside" example of a matrix for which eigenvectors are "relatively simple," but which also has a "relatively simple" vector for which the magnitude is preserved, although not necessarily the direction. Your reply at the end of the program was very insightful, and I think it opens an interesting are of exploration.

Juan Tolosa Richard Stockton College, NJ

It looks like the Online Help has now been updated for Maple 2015.

Here's the What's New in Maple 2015 Overview, which is similar to the product pages currently available here.

But we can now also link to and view the online versions of full help pages of new or enhanced commands or example worksheets. For example, dataplot.

acer

@Markiyan Hirnyk   I try not to use this package, as I think the results are not reliable enough. Here is the example, where instead of the three real roots it finds only one, despite the hint to look for the three roots:

restart;

DirectSearch:-SolveEquations(100^x=x^100, AllSolutions, solutions=3);

There are many other examples, particularly in discrete optimization in which it returns false results. Here is one example (well-known to you).

I would like to pay attention to an article " Sums and Integrals: The Swiss analysis knife " by Bill Casselman, where the Euler-Maclaurin formula is discussed.  It should be noted all that matter is implemented in Maple through the commands bernoulli and eulermac. For example,

bernoulli(666);

eulermac(1/x, x);

,

eulermac(sin(x), x, 6);

BTW, as far as I know it, this boring stuff is substantially used in modern physics. The one is considered in

Ronald Graham, Donald E. Knuth, and Oren Patashnik, Concrete mathematics, Addison-Wesley, 1989.

The last chapter is concerned with the Euler-MacLaurin formula.


           

Here are some tips and tricks that will help you get the most out of Maple 2015, covering from short cuts to how to use the newest features.

1. Whenever you are asking yourself “..but how do I do it?”, just type ?Portal+Enter, and you will access the Maple Portal, which will give you a complete guide on how to do things.
   
   
2. If you want to implement 1 of the 300 tasks that Maple offers in a syntax-free way, like Completing the Square, just follow this path: Tools≻Tasks≻Browse.
   
   
3. Type Ctrl+F2 or Command+F2 and the Quick Reference window with shortcut keys and other information about working with the Maple interface will pop up.
   
   
4. If you need quick help with a specific mathematical function, click or highlight the function + F2 and a Help box that contains a summary of the basic characteristics of the function will pop up.
   
   
5. If you have installed the Excel Add-in and you want to perform some Maple commands within Excel, make sure to enable the Maple add in by following this path: Excel’s Tool Menu>Add-Ins>Select Maple Excel Add-in Box> OK
   
   
6. Export Maple’s data into Excel by right clicking and choosing ‘Export As’>Excel.
   
   
7. Instead of having to copy-paste your Maple information into a Power Point Presentation, just turn the slideshow mode on by pressing F11. This way you will have an interactive presentation that holds all the live plots and embedded components that Maple offers.
   
   
8. Whenever you want to create interactive mini-applications that can be used to explore the parameters of any arbitrary Maple expression, such as a plot, mathematical equation, or command use the Exploration Assistant. Do this by either right-clicking +Explore from the context-sensitive menus, or by calling the Explore command.
   
   
9. Save time while computing mathematical expressions by calling the equation label instead of having to re-type the equation. Do this by pressing CTRL+L and then input the number that identifies the equation.
   
   
10. Reference mathematical equations or expressions from other documents. First, determine which label is associated with the equation you want. In the main document, select "Insert" > "Reference". From the file dialog, select the file containing the expression. Then select the equation reference number of your equation from the list that appears.
    
    
11. In Maple, the letter "e" entered using the keyboard does not represent the exponential function. The exponential function can be entered using command completion (Ctrl+Space or ESC) or the "exp(a)" item in the Expression Palette (Standard interface only). The exponential can also be entered as:          

    > exp(x)
    
    

12. With Maple 2015 you can now access data sets from various built-in and online data sources. This package is able to access time series data from the data aggregator Quandl, as well as locally installed data from countries and cities. To learn more, click here.
    
    
13. Whenever you assign plots to a variable name, p:=[plot(sin(x)), plot(cos(x))] a thumbnail of the plot will appear instead of the code.
    
    
14. Save time when inputting existing or personalized units. Just click CTRL+SHIFT+U and type the desired units you want.
    
    
15. With Maple 2015 you can now zoom in or out just by pressing CTRL+SCROLL or CTRL+ place two fingers on the pad and move them up to zoom in or down to zoom out.
    
    
16. Convert a Maple Worksheet into Microsoft Word: This can be done using the Export to HTML feature.
    
    1. Prepare your worksheet as you would like it to appear in the document.
    2. From the "File" menu in Maple, select "Export As ..." > "HTML".
    3. Give the HTML file a name, "output.html" for example.
    4. When the export has completed, start Word, and open the HTML file. If you used "output.html" as the name to save the file as, open the file called "output1.html" into Word.
    5. From the "File" menu in Word, select "Save as Word Document" to save the file. You now have a Word document which contains the content of your Maple worksheet.
       
       Note: this procedure will work with any Word Processing program that can open an HTML document.
       
       
17. Change Maple’s default input from 2D to 1D:
    1. Open the Tools > Options... menu (Maple > Preferences on a MACINTOSH machine).
    2. Select the Display tab
    3. From "Input Display" menu select Maple Notation
    4. Press the Apply to Session button to make the change take effect for the current Maple session.
    5. Press Apply Globally to have the change take effect permanently. Maple will need to be restarted if you choose Apply Globally for the changes to take effect.
       
       You may download a set of instruction on how to change your 2D interface to the “Classic” Style here: ftp://public.maplesoft.com/miscellaneous/ChangeToClassicInterface.pdf

We hope that you find this list helpful. Please feel free to add any of your tips or techniques to this post, or to create your own new topic.

Last month, we received a very kind note from a recipient of one of our sponsorships. Maplesoft sponsors several academic and commercial events throughout the year, providing free copies of Maple or MapleSim to lucky attendees. Audrey was one of the winners of the Elgin Community College Calculus Contest, where she won a copy of Maple. Here’s what she had to say:

Thank you so much for the Maple license.  I have become familiar with Maple during the last school year.  At first the commands were like Chinese to me and I had a rough time getting anything done, but once I made a connection between the commands and what they were doing it was a lot easier.  Even without former knowledge of computer programing, the commands are increasingly intuitive.  Maple has been a huge help to me doing my homework and projects, and even as I was studying for the competition it was useful for checking my answers.  Another reason that I love Maple is that it provides visuals for the difficult concepts we learned in class, such as shell method in Calc II and mixed partial derivatives in Calc III.  I enjoy math, but I thank that Maple has enriched my experience along the way.

Thank you again for your generous gift, 

~Audrey~

It’s always nice to hear how students and professionals alike are succeeding with the help of Maple. If you’d like to share your experience, please send an email to customerservice@maplesoft.com or post it here on MaplePrimes.

You are teaching linear algebra, or perhaps a differential equations course that contains a unit on first-order linear systems. You need to get across the notion of the eigenpair of a matrix, that is, eigenvalues and eigenvectors, what they mean, how are they found, for what are they useful.

Of course, Maple's Context Menu can, with a click or two of the mouse, return both eigenvalues and eigenvectors. But that does not satisfy the needs of the student: an answer has been given but nothing has been learned. So, of what use is Maple in this pedagogical task? How can Maple enhance the lessons devoted to finding and using eigenpairs of a matrix?

In this webinar I am going to demonstrate how Maple can be used to get across the concept of the eigenpair, to show its meaning, to relate this concept to the by-hand algorithms taught in textbooks.

Ah, but it's not enough just to do the calculations - they also have to be easy to implement so that the implementation does not cloud the pedagogic goal. So, an essential element of this webinar will be its Clickable format, free of the encumbrance of commands and their related syntax. I'll use a syntax-free paradigm to keep the technology as simple as possible while achieving the didactic goal.

Notes added on July 7, 2015:

• This live webinar has already been presented and you can view the recording here: http://www.maplesoft.com/webinars/recorded/featured.aspx?id=922
  
  
• For one reaction to a point made in the webinar, see this post.

We find recent applications of the components applied to the linear momentum, circular equations applied to engineering. Just simply replace the vector or scalar fields to thereby reasoning and use the right button.

Momento_Lineal_y_Circular.mw

(in spanish)

Atte.

L.AraujoC.

A small piece of code for fun before the weekend, inspired by some recent posts:

http://www.johndcook.com/blog/2015/06/03/mystery-curve/

http://mathlesstraveled.com/2015/06/04/random-cyclic-curves-5/

http://www.walkingrandomly.com/?p=5677

cycler := proc(t, k, p, m, n, T) local expr, u, v;
  u := exp(k*I*t);
  expr := exp(I*t) * (1 - u/m + I*(u^(-p))/n);
  v := 1 + abs(1/n) + abs(1/m);
  plots:-complexplot( expr, t = 0 .. T, axes = none,
                      view = [-v .. v, -v .. v] );
end proc:

cycler(t, 5, 3, 2, 3, 2*Pi);

And that can be made into a small application (needs Maple 18 or 2015),

Explore( cycler(t, k, p, m, n, T),
         parameters = [ k = -10 .. 10, p = -10 .. 10,
                        m = -10.0 .. 10.0, n = -10.0 .. 10.0,
                        [ T = 0 .. 2*Pi, animate ] ],
         initialvalues = [ k = 5, p = 3, m = 2, n = 3, T = 2*Pi ],
         placement = left, animate = false, numframes = 100 );

cyclerapp.mw

The animation of parameter T from 0 to 2*Pi can also be fun if the view option is removed/disabled in the call to complexplot. (That could even be made a choice, by adding an additional argument for calling procedure cycler and an accompanying checkbox parameter on the exploration.)

acer

Maplesoft will be hosting the 2015 Maple T.A. User Summit this June 15 - 17 in New York City. Don’t miss this opportunity to learn about new trends in online education while networking and socializing with fellow educators and Maple T.A. users in the city that never sleeps!

We are happy to announce that the schedule has been finalized! The event will include keynote presentations, talks and discussions from users and Maplesoft staff, training sessions, a welcome reception and a boat cruise around New York City, 

If you'd like ot sign-up but still haven't - don't hesitate to do so today using the following link: https://webstore.maplesoft.com/taconference/register.aspx .

I hope to see you there!

Jonny
Maplesoft Product Manager, Maple T.A.

Developed and then implemented with open code components. It is very important to note this post is held for students of civil engineering and mechanics. Using advanced mathematical concepts to concepts in engineering.

Metodos_Energeticos_full.mw

(in spanish)

Atte.

L.Araujo.C

On a related post here, I have posted some screen shots for fun. Today is a special day for me in particular, because I passed my PhD viva with minor corrections! If you look colsely on the date, by coincidence, it was TWO years ago exactly. So I feel like I maybe 'update' the old post with a bit more, to make it V2.

So here comes the new screent shots: (Bear with me for now, I will post some actual content about how I feel when using Maple later.)

Alright, enough with the screent shots. If you click to see the big picture, you would notice that I have managed to get a complete sets of Maple versions from Maple V to Maple 2015!

===============================================

I am originally from Shanghai and I have always done well in Maths and science. I first heard about Maple in 2004 from a Maths teacher. He introduced me to this software. I played just a couple of times with Maple 6 on his computer, to get a first impression. At that time, all I could say was, hmmm insteresting.

In 2006, I went to UK for a foundation course. That was the first time I was actually taught how to use Maple 9.5. So I had access to it on the university computers. I discovered a lot about Maple and used it for ALL my Maths homeworks. Yes, I am lazy. I hope that Maple can do everything!

Then I got fascinated with Maple. Unfortunately, Maple was not taught in the undergraduate course. But there are materials for self-learning. By that time, I had become quite good with most of the contents in those materials. So I went look for other things to try. I went to ask my foundation course classmates, to see if they have anything on Maple and if they needed the help. One of them was at Imperial College and the questions were a bit chanlleging (finally!). That was the time I first met MaplePrimes! Hello, how have you been?

After that I have never stopped using Maple or Maple Primes. I may have been quite on the forum, but I was there. My PhD research rely hugely on the ability to compute the symbolic rank of certain matrices as well as the ability to simplifiying complicated expressions using siderules.

I do not want to talk too much about the technical details when using Maple. But I do

THANK Maple and it creators, developers and all relavent staff around the globe.

THANK Maple Primes and all its users. Some users have been particulary helpful!

Lastly, I dont have any other words to say other than this:

It's been really fun and enjoyable. Thanks!

This post is related to the this thread

The recursive procedure  PosIntSolve  finds the number of non-negative or positive solutions of any linear Diophantine equation

a1*x1+a2*x2+ ... +aN*xN = n  with positive coefficients a1, a2, ... aN .  

Formal parameters: L is the list of coefficients of the left part, n  is the right part,  s (optional) is nonneg (by default) for nonnegint solutions and  pos  for positive solutions.

The basic ideas:

1) If you make shifts of the all unknowns by the formulas  x1'=x1-1,  x2'=x2-1, ... , xN'=xN-1  then  the number of positive solutions of the first equation equals the number of non-negative solutions of the second equation.

2) The recurrence formula (penultimate line of the procedure) can easily be proved by induction.

The code of the procedure:

restart;

PosIntSolve:=proc(L::list(posint), n::nonnegint, s::symbol:=nonneg)

local N, n0;

option remember;

if s=pos then n0:=n-`+`(op(L)) else n0:=n fi;

N:=nops(L);

if N=1 then if irem(n0,L[1])=0 then return 1 else return 0 fi; fi;

add(PosIntSolve(subsop(1=NULL,L),n0-k*L[1]), k=0..floor(n0/L[1]));

end proc:

Examples of use.

Finding of the all positive solutions of equation 30*a+75*b+110*c+85*d+255*e+160*f+15*g+12*h+120*i=8000:

st:=time():

PosIntSolve([30,75,110,85,255,160,15,12,120], 8000, pos);

time()-st;

                                       13971409380

                                             2.125

To test the procedure, solve (separately for non-negative and positive solutions) the simple equation  2*x1+7*x2+3*x3=2000  in two ways (by the  procedure and brute force method):

ts:=time():

PosIntSolve([2,7,3], 2000);

PosIntSolve([2,7,3], 2000, pos);

time()-ts;

                47905

                47334

                 0.281

ts:=time():

k:=0:

for x from 0 to 2000/2 do

for y from 0 to floor((2000-2*x)/7) do

for z from 0 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k; 

k:=0:

for x from 1 to 2000/2 do

for y from 1 to floor((2000-2*x)/7) do

for z from 1 to floor((2000-2*x-7*y)/3) do

if 2*x+7*y+3*z=2000 then k:=k+1 fi;

od: od: od:

k;

time()-ts; 

                   47905

                   47334

                   50.063

Another example - the solution of the famous problem: how many ways can be exchanged $ 1 using the coins of smaller denomination.

PosIntSolve([1,5,10,25,50],100);

                        292

 Number-of-solutions.mw

 Edit.  The code has been slightly edited