Category: Differential Equations

Eulix: A Stiff ODE Solver in Maxima

Look what I found in the share/contrib directory:

There isn’t any documentation available from texinfo:

Browsing around on the web I found this documentation page.

There are also many examples given in the Eulix_*.mac files in share/contrib/Eulix . Here’s what is says at the top of Eulix.mac:

* (C) Helmut Jarausch
 *
 * @version 1.0 Date August 29th, 2017
 *
 * @author Helmut Jarausch
 *         RWTH Aachen University

 * @contents
 *  Eulix is an linearly implicit (by default) extrapolated Euler method.
 *  It uses variable step sizes, variable orders and offers dense output.
 *  Eulix can solve systems like  M y'(t)= f(t,y)  where M is a constant mass
 *  matrix. This enables the solution of differential algebraic equations (DAE)
 *  of index 1, as well.
 *
 *  I have made a simplified version of  SEULEX  published in
 *
 *  E. Hairer, G.Wanner
 *  Solving Ordinary Differential Equations II
 *  Stiff and Differential-Algebraic Problems
 *  Springer 1991

A Test Problem: The Ball of Flame

To test drive the method, I’ll retrace Cleve Moler’s steps in this excellent piece about numerical solvers for stiff equations.

We’ll consider the “ball of flame” equation of Lawrence Shampine:

\dot{y}=y^2-y^3, \;\;\;\;\;\; y(0)=\delta

on the time interval 0 \le t \le 2/\delta .

The equation is meant to be a simple model of flame propagation, but it could also describe what happens to non-stiff solvers for small values of \delta : they go down in a ball of flames!

Notice that the basic calling sequence for Eulix() is the same as for rkf45(). The first thing to observe is that when applied to the ball of flame, Eulix() uses less time and less steps than rkf45():

Plots of the solution points show that Eulix() makes this look like an easy problem, using very few points for the horizontal parts of the solution, and concentrating more points at the sharp turns:

The non-stiff variable step Runge-Kutta-Fehlberg method rkf45() makes this look hard:

Zooming in shows how hard rfk45() is working on the horizontal part of the solution. It is doing its job to maintain the error within the default tolerance 1e-6, but it’s using a lot of steps and function evaluations to do that:

And one final thing to know about Eulix(): The output can appear very sparse, as we see in the region t> 10000 in the figure above. There’s an easy option to specify dense output (likely computed using splines and the existing solution points). Add the desired output interval as a fourth entry in the time specification list: [t,tstart,tfinal,output_interval]

Thank you Helmut Jarausch!

Quotes, global variables, and ODE solvers.

Computed solution: room temperature in a model of a thermostatically switched air conditioning system

Every winter term in my introductory ordinary differential equations class, I ask each of my students to choose from the scientific literature a paper that puts ODE models at the center of the work, and to reproduce the author’s results. For the past 5 years we’ve used Maxima for that purpose, and throughout the course for both symbolic and numerical tasks.

A perennial favorite is “Temperature Models for Ware Hall” by J. K. Denny and C. A. Yackel, The College Mathematics Journal, Vol. 35, No. 3 (May, 2004), pp. 162-170 http://www.jstor.org/stable/4146891 .

We have a first order ODE that combines two Newton’s law of cooling processes: the room is in thermal contact with the hot outside air at 90F and with cooled air at 60F coming from the air conditioner.

\frac{dT}{dt}=0.2(90-T)+0.2S(T)(60-T)

The idea is to model an air conditioning system that is controlled by a thermostatic switch. In the model above the switch function S(T)=1 when the cooling sytem is running, becomes zero if the system is running and the temperature reaches 67F, becomes 1 if the system is not running and the temperature reaches 73F. The behavior of the room temperature is shown in the figure at the top of this piece.

There might be another way to model this (I’d be delighted to learn what that is!) but here we treat the thermostat with a global variable that can be read and changed within a numerical solver as part of the right-hand-side function.

The question in this post is: how to implement a global flag like that in Maxima? Below are four variants involving various combinations of the Maxima quote operator. The first two give the hoped-for behavior and were used to produce the figure above. The last two don’t seem to have access to the global flag and the cooling system runs without ceasing. I include a graph of that incorrect computed solution at the bottom of this post.

The Details

  • In the first case, the switch function is included in the right-hand-side function with a single quote.
  • In the second, the switch in included with quote-quote, but the order in which the right hand side and switch is important: the switch function must be defined AFTER the right hand side function in which it is called
  • In the third, the switch in included with quote-quote, but the order from above is reversed and the global flag seems undetected.
  • In the fourth, the switch is included without quote or quote-quote and the behavior is same as quote-quote in the second and third cases.

I don’t know why each works or doesn’t! I’m hoping you do and can tell me 🙂 Thanks in advance.

/*This works as intended with single-quote on ThermSwitch
    ThermSwitch can be defined either before or after ff */
kill(all);
globalAC:1;

ThermSwitch(T):=block(
	      if T<67 then globalAC:0
	         else if T>73 then globalAC:1,
	      return(globalAC))$
		  
ff:.05*(90-T) + .2*'ThermSwitch(T)*(60-T);

s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));
/*This works as intended with double-quote on ThermSwitch
    ThermSwitch  defined  AFTER ff */
kill(all);
globalAC:1;

ff:.05*(90-T) + .2*''ThermSwitch(T)*(60-T);

ThermSwitch(T):=block(
	      if T<67 then globalAC:0
	         else if T>73 then globalAC:1,
	      return(globalAC))$
		  
s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));
/*This does NOT work as intended with double-quote on ThermSwitch
    ThermSwitch  defined  BEFORE ff */
kill(all);
globalAC:1;

ThermSwitch(T):=block(
	      if T<67 then globalAC:0
	         else if T>73 then globalAC:1,
	      return(globalAC))$
		  
ff:.05*(90-T) + .2*''ThermSwitch(T)*(60-T);

s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));
/*This does NOT work as intended with unquoted ThermSwitch before ff
   it does work as intended with unquoted ThermSwitch after ff  */
kill(all);
globalAC:1;

ThermSwitch(T):=block(
	      if T<67 then globalAC:0
	         else if T>73 then globalAC:1,
	      return(globalAC))$
		  
ff:.05*(90-T) + .2*ThermSwitch(T)*(60-T);

s:rk(ff,T,78,[t,0,50,.1])$
wxdraw2d(points(makelist([p[1],p[2]],p,s)));
Incorrect Computed solution: room temperature in a model of a thermostatically switched air conditioning system — the thermostatic switch in the model isn’t functioning correctly.

Multiple Plots in an Animated GIF in Maxima

Here it is: Roughly the animation I wanted to create. A harmonic oscillator, with the corresponding trajectory in phase space and a bead-on-wire projection of the trajectory onto the potential energy curve.

I wanted to insert this in a powerpoint slide I would be using for an online lecture, so an animated GIF seemed like a good option. I started with an example from Dirk Mittler’s blog. (link no longer working)

What I really wanted to do was include several axes like in these examples from my earlier post at TheMaximalist, which created multiple scenes using several calls to gr2d() and then drew the scenes in a single call to draw(). The trouble I encountered was that the collection of time-synchronized multiple scenes weren’t rendered simultaneously in the animated GIF — the animation flickered back and forth between the scenes.

So I turned off the axes and made my own triplet of axes in a single scene using parametric lines and including the needed vertical offsets.

One key thing that took me a while to find: In my Windows installation (Maxima 5.44 and wxMaxima 20.06.6), no matter what I specified for a path in the file_name argument to draw(), the animated GIF file ended up getting written in C:\maxima-5.44.0

I’ve pasted the code below.

scenes: []$

for i:0 thru 61 do (
    scenes: append(scenes,
            [gr2d(proportional_axes=xy,
    noframe,
    ytics=false,xtics=false,nticks=100,
    /* oscillator 2 units below */
    color=black, 
      parametric(t,-2,t,-1.2,1.2),
    /* phase space 2 units above */
    color=black, 
      parametric(t,0,t,-1.2,1.2),
      parametric(0,t,t,-1.2,1.2),
    color=blue,
    parametric(cos(t),sin(t),t,0,2*%pi),
    /* potential */
     color=black, 
      parametric(t,2,t,-1.2,1.2),
      parametric(0,2+t,t,-.1,1),
    color=blue,
     explicit(2+x^2,x,-1.2,1.2),
    color=red,point_type=filled_circle,
        points([cos(i/10),cos(i/10),cos(i/10)],[-2,-sin(i/10),2+cos(i/10)^2]))]
        )
)$

draw(
    delay = 10,
    file_name = "oscillator_phasespace_potential",
    terminal = 'animated_gif,
    scenes
)$

An improved Maxima function for inverse Laplace transform

Maxima has a fairly serviceable Laplace transform utility built-in.  Here’s an example from the popular ordinary differential equations book by Blanchard, Devaney and Hall:

laplace1

Trouble arises when we look at discontinuous forcing functions, which is especially sad because it seems to me that’s what makes it worthwhile to spend time with Laplace transforms.  In particular, the inverse transform function ilt() fails on the Heaviside Function and Dirac Delta, even though the built-in laplace() treats them correctly:

laplace2

So, I’ve written an alternative inverse Laplace function laplaceInv() that fixes that problem:

laplace3

Here are a few differential equation solutions to show how the new function behaves:

A second-order linear equation with a constant forcing function that vanishes at t=7

laplace6

laplace7

A second-order equation with two impulsive forces:

laplace4

laplace5

The new Laplace transform Maxima function can be downloaded here.

An undocumented synonym for diff() in Maxima

Today, a student turned in some Maxima work for my class. I discovered he had successfully used the command derivative() in place of diff() with seemingly identical results.  I verified that the same thing works in several versions of Maxima I have installed on my windows computer.  Who knew?

diffderivative

logabs, logarc: How to make integrate() return what you expect

logarc

In the back of my calculus book there is a table of famous integrals.  Here’s integral number 21 in that table:

20171004_113109

From Maxima integrate(), I get

integral21Maxima

 

What’s going on?

Both forms give a workable antiderivative for the original integrand:

integral21Maxima2

Furthermore, we believe that both forms are correct because of this helpful identity for hyperbolic sine:

 {\rm asinh}(z)=\ln(z+\sqrt{1+z^2}).  

Turns out (thanks to a Barton Willis for pointing me in the right direction) there’s a variable logarc that we can set to make Maxima return the logarithmic form instead of hyperbolic sine:

integral21Maxima3

I haven’t yet encountered cases where this would be a bad idea in general, but I’ll update this if I do.

logabs

In the first week of my differential equations course, we study methods of direct integration and separation of variables.  I like to emphasize that the absolute values can lend an extra degree of generality to solutions with antiderivatives of the form

\int \frac{1}{u}\;du = \ln |u|+ C.   

As an example, for the initial value problem

 y' = \frac{x}{1-x^2}    ,  y(0)=1    ,

it is convenient for treating all possible initial conditions (x \ne \pm 1) in one step to use the antiderivative

y = -\frac{1}{2} \ln | 1-x^2 | + C.   

However, Maxima omits the absolute values.

logabs1.PNG

For this case, we could consider only the needed interval -1<x<1, but still…

Turns out we can set the Maxima variable logabs to make integrate() include absolute values in cases like this:

logabs2.PNG

But then later in the course, I saw that logabs also impacts the Ordinary Differential Equation solver ode2().  I encountered an example for which Maxima, in particular solve() applied to expressions involving absolute value,  didn’t do what I wanted with logabs:true

For the logistic equation

\frac{dP}{dt} = kP\left ( 1-\frac{P}{P_c} \right )    ,  P(0)=P_0   

we expect that by separating variables we can obtain the solution

P(t) = \frac{P_0 P_c}{(P_c-P_0)e^{-kt}+ P_0}.   

Here’s what happens when we use ode2() with and without logabs:true:

logistice2

logistic3

 

 

Student Projects in Maxima Vol 1, Num 4: ODE Systems for Epidemiological Models

Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences.

Volume 1, Number 4 (2017) is devoted to Systems of Ordinary Differential Equations for SIR models in epidemiology.

 

Contents:

Emenheiser, Anna, Virus Dynamics and Drug Therapy

Radermacher, Erin, A Model of the 2014 Ebola Outbreak

Rafai, Sagar, Analysis of Communicable Disease Models

 

Student Projects in Maxima Vol 1, Num 3: ODE Systems for Interacting Population Models

Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences.

Volume 1, Number 3 (2017) is devoted to Systems of Ordinary Differential Equations for non-chaotic predator-prey and other interacting population models

 

Contents:

Bhullar, Abhjeet,  Lions, Wildebeest and Zebras

Goodin Aunic, Parasitic Nematodes in Grouse

DeFore, Brandon, Breeding Suppression in Predator-Prey

Jerez, Emilio, Predator-Prey with Hawk and Dove Tactics

Bontrager, Eric, Predator-Prey with Prey Switching

Beatty, Ethan, Analysis of Logistic Growth Models

Rice, Gabriel, Pharmacokinetic Mechanism of Ethanol-Benzodiazepine Interactions

Kay, Ian, Radioactive Isotopes

Wile, Jessica, Ebola in Western Lowlands Gorillas

Bailey, John,  Kinetic Modeling for Interconnected Reactions

Piet, Joe, Elephant and Tree Population Dynamics

Kim, Judy, A Model for West Nile Virus

Kamalaldin, Kamalaldin, Infected Prey in Polluted Environment

Park, Kayla, Ratio-Dependent Predator Prey Model

Lundy, Liam, Predator Prey in a Single Species with Cannibalism

Orwin, Michael, Interactions of Model Neurons

Schultz, Pete, Photosynthetic Oscillations

Wadhwa, Raoul, Dynamics in an Inflammation Model

Del Olmo, Ricardo, Population Growth and Technological Change

Network Model of Cocaine Traffic in Spain Edited

Kill, Sean, A Logistic Model with Varying Carrying Capacity 

McFadden-Keesling, Sophia,  Predator Prey Model with a Refuge

Samson, Tanush,  Lions and Wildebeest Model

Morales, Zach,  Game Theory Models

Student Projects in Maxima Vol 1, Num 1: Systems of ODEs and Chaos

Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences.

Volume 1, Number 1 (2017) is devoted to Chaotic Systems of Ordinary Differential Equations.

Contents:

Thornburg, Eric, Simple Pendulum and Chaos

Bhimani, Kevin, Chaos in a 3-Species Food Chain

York, Lily, A Lattice Model of Epilepsy

Andrews, Steven, 3D Chaotic Model

Rutledge, Tim,  A Model of Neuronal Bursting

 

 

Student Projects in Maxima Vol 1, Num 2: Systems of ODEs with Impulses and Switching Functions

Student Projects in Maxima is an online project to disseminate work by undergraduate students using Maxima to reproduce published results in the sciences and social sciences.

Volume 1, Number 2 (2017) is devoted to Systems of Ordinary Differential Equations with time-dependent forcing terms, non-continuous inputs, and forcing terms that require knowledge of a past state of the system.

Contents:

Chumley, Qynce, An Office Heating and Cooling Model

Williams, Nick, Baseball Pitch Dynamics

Rizzolo, Skylar, BAC Model for alcohol consumption

Barth, Eric, A model of bladder bacteria proliferation in prostate disease

Barth, Eric, A Model of Pulse Vaccination Strategyi