# Evaluation problems with NDSolve using the fixed step ExplicitEuler-Method

## Motivation

I am using Mathematica to implement System Dynamics (SD) models, e.g. differential equation models as applied to management/economic problems following the "methodology" of Jay W. Forrester. Widespread specialist software like Ventana's Vensim use the simple explicit Euler method as the advantages in ease of modeling/programming are assumed to outweigh the numerical disadvantages of this, rather simple, fixed step solution method.

I am following that trail here -- also for other reasons than pure conformity. In other words, let us simply assume that for some reasons the explicit Euler method is to be applied to solve the ODEs.

## Problem

In a simple model I am modling the payment to an account as a pulse process so that the development in the stock variable (eg. the account) should follow a step-function. Naturally in continuous time a sequence of DiracDelta-Functions would be the way to go but I have found them to be incompatible with the explicit Euler method so far.

## Model

Needs["DifferentialEquationsNDSolveUtilities"];

pulseTrain::usage = "pulseTrain[start,width,tbetween,end] will return a
a function of time which will provide a pulse-sequence starting at start
and repeating in intervals of length tbetween until end-time. Each pulse
will have the length indicated by width. The amplitude of the pulse
is 1.";
pulseTrain[ start_, width_, tbetween_, end_] := If[
tbetween <= width,
(* then *) Function[time, Piecewise[{{1, start <= time <= end}}, 0]],
(* else *) Function[time, Piecewise[
Table[
{ 1, ti <= time < ti + width},
{ti, start, end, tbetween}
],
0
]
] (* endif *)
]

b[t_] := pulseTrain[ 0 , 1/32 , 1 , 10 ][t]; (* raw pulse *)

g[t_] := 32 b[t] (* amplitude to be 1/(width of pulse) *)

sol = NDSolve[
{
(* Net Flows *)
x1'[t] == g[t],
x2'[t] == b[t],
(* Initial Stocks *)
x1[0] == 100,
x2[0] == 100
},
(* stocks *)
{ x1, x2 },
(* time range *)
{t, 0, 10},
Method -> "ExplicitEuler",
StartingStepSize -> 1/32
];


## Analysis

Here is the plot for the functions g(t) and b(t) defining the pulse processes, where g(t) should integrate to +10 over the interval [0,10].

Using EulerIntegration as indicated in the NDSolve-Function the correct flow g(t) will integrate to 320 giving an end value for the account of 420. The raw pulse b(t) will suprisingly give the result of 110 which one should expect for the integration of g(t):

The correct result for x1(t) (e.g the integration of g(t)) can be obtained using StartingStepSize-> 1/33 (* or smaller *) or by using Method -> "ExplicitRungeKutta" (* without any StartingStepSize*).Note, that NDSolve will not find the correct solution without any options given.

I also noted that, although I had used a fixed step method (eg. Euler-Integration), StepDataPlot[ sol, PlotRange -> Full, PlotTheme -> "Detailed"] reveals the following plot, which surprises me, as it reveals alternating step sizes :

What is going on here: Why is the simple Euler method with a sufficient resolution (eg. StepSize == Width of Pulse) not showing the correct result and why is this fixed step method using a variable time step?

Remark: Needless to say that the equivalent model in Vensim using identical parametrization (eg. pulseWidth = StepSize = 1/32 ) returns the correct end value of 110.

• I would dare say that "the advantages in ease of modeling/programming" are totally outweighed by using NDSolve. All the work of programming many methods has been done for you. I believe this remark is addressed to those who would have to write their own program and have limited experience with programming and numerics (or have limited time to develop a sophisticated approach). Therefore, I would suggest that you let NDSolve choose its own method, unless after doing so, it seems to give bogus results. – Michael E2 May 18 '15 at 18:33
• @Michael E2: So I had thought also but if you let NDSolve choose its own method (eg. comment out StartingStepSize and Method options) it will not give the correct result which I do find disappointing somehow. – gwr May 18 '15 at 18:37
• @Michael E2: I am sticking to EulerIntegration with fixed step sizes because I am employing the Bayesian Particle Filter to identify parameters in my models and -- so far ;-) -- the resampling method build upon a fixed step markov chain (no continuous particle filter yet). I would really "love" to forget all about StepSize and simply "build in continous time" like a "purist"... – gwr May 18 '15 at 18:40
• try defining your b and g as follows: b[t_?NumericQ]:=pulseTrain[0,1/32,1,10][t];g[t_?NumericQ]:=32 b[t] – chuy May 19 '15 at 14:32
• @gwr My (last) solution was to use "DiscontinuityProcessing" -> False with a step size of 1/32, which is what happens when you have function protected by NumericQ. The DiracDelta was a gratuitous addendum, which cannot work with regular Euler scheme, without approximating it by a finite pulse, as you do with your b[t]. I wanted to show that discrete events can be modeled and one does not have to stick Euler schemes. By "liking" I just meant you just ignored the primary solution and criticized the Dirac delta one. – Michael E2 May 19 '15 at 15:31

## 2 Answers

First, you're not using a fixed step method. (An Euler scheme may be applied to any step size and to one that varies.) To get a true fixed step method you have to turn off "DiscontinuityProcessing" when you have a discontinuous ODE; otherwise, NDSolve will try to adapt the steps to account for the discontinuity. The "DiscontinuityProcessing" stage resets the step size when a discontinuity is detected. Turning it off gives a straightforward, fixed-step Euler scheme. Second, it is not clear why NDSolve seems to use the value 1 or 32 for the value of b[t] or g[t] respectively in your setup. It certainly seems like a bug.

{sol} = NDSolve[
{x1'[t] == g[t], x2'[t] == b[t], x1[0] == 100., x2[0] == 100.},
stocks = {x1, x2},
{t, 0, 10},
Method -> {"FixedStep", Method -> "ExplicitEuler",
"DiscontinuityProcessing" -> False}, StartingStepSize -> 1/32];

Plot @@ {Through[stocks[t]] /. sol,
Flatten[{t, First[stocks]["Domain"] /. sol}]}


Addenda

While the Euler method is a bit of a dinosaur and not very accurate, you may feel that you cannot abandon it. Nevertheless, there are other ways to model a discontinuous process. Obviously an Euler method cannot by itself handle a discontinuity like a Dirac delta-function. One would have to approximate it with a unit-area box function, which is basically what the OP has done. To take advantage of the sophisticated methods in NDSolve, you cannot hobble it with the Euler method.

Here is an example of using DiracDelta.

{sol} = NDSolve[
{x'[t] == Sum[DiracDelta[t - t0], {t0, 0, 10}], x[0] == 100.},
stocks = {x}, {t, 0, 10}];

Plot @@ {Through[stocks[t]] /. sol,
Flatten[{t, First[stocks]["Domain"] /. sol}]}


Another method using events:

{sol} = NDSolve[{x1'[t] == 0, x1[0] == 100,
WhenEvent[Mod[t, 1] == 0, x1[t] -> 1 + x1[t]]},
stocks = {x1}, {t, 0, 10}, StartingStepSize -> 1/32];


Both produce a nicer graph:

• Thanks for showing the DiracDelta- function but I had mentioned and excluded that possibility in my post: Naturally in continuous time a sequence of DiracDelta-Functions would be the way to go but I have found them to be incompatible with EulerIntegration method so far. – gwr May 18 '15 at 19:32
• @gwr I've undeleted my answer. I added a little explanation about the discontinuity processing and about Dirac, and one other method for handling discrete events. If you are stuck with Euler, they won't be of use to you, but I think others who visit the site might find the alternatives helpful. – Michael E2 May 19 '15 at 15:35
• Thank you. I corrected your WhenEvent solution to show x1'[t] == 0. Also note, that one can very well use this kind of event management with Euler - in fact I am using it to do the sampling-importance resampling with the particle filter. The latter is the reason I (believe) to need a fixed time step as equation and measurement noise are defined as discrete stochastic processes. Maybe I will not need fixed steps after all? – gwr May 19 '15 at 16:39
• @gwr The error was from playing with your system and copying the wrong code -- thanks for fixing it. Certainly discrete processes can (or ought to) be modeled with fixed steps. But NDSolve can handle discrete events, too. I suppose it's not clear that there is a choice that works better in all cases, especially given the difficulty you encountered. – Michael E2 May 19 '15 at 16:57
• FWIW, some Q&A with a manual implementation of Euler's method: mathematica.stackexchange.com/q/22042, mathematica.stackexchange.com/q/41783, mathematica.stackexchange.com/q/22413, mathematica.stackexchange.com/a/30577 – Michael E2 May 19 '15 at 17:01

## Using Event-Handling and Euler method

Expanding the solution given by Michael E2 one might figure that using EventSeries is coming closest to the "real thing" which after all is a series of discrete events:

eventTimes = Range[ 0, 10, 1 ]; (* example *)

isEventTimeQ[ t_?NumericQ ] := Piecewise[
Table[ { 1 , t == eventTime }, {eventTime, eventTimes}],
0
]

cashflow = EventSeries[
RandomInteger[{-10, 20}, Length @ eventTimes],
{eventTimes}
]  (* example *)

{sol} = NDSolve[
{
(* net flows *)
account'[t] == 0, (* or later on: compoundingRate x account[t] *)
(* initial stock values *)
account[0] == 100,
(* discrete events *)
WhenEvent[ isEventTimeQ[t] == 1,
account[t] -> cashflow[t] + account[t]
]
},
(* stocks *)
account,
(* time range *)
{t, 0, 10},
(* options *)
StartingStepSize -> 1/4, (* event times must coincide with eval-times! *)
Method -> "ExplicitEuler"
];


## Results

Thus with a series of discrete cash flows:

On gets a true step function that might later on be mixed with continous compounding in an expanded model:

As far as I can see this solution only works with the Euler method as the isEventQ[t] function only works inside the WhenEvent framework if the event times match evaluation time steps. But of course one might then use the modulus-formulation inside WhenEvent, as Michael E2 has done, and use MissingDataMethod -> {"Constant",0} to set up an irregular EventSeries to match the regular sampling interval in that case.

• If it helps, this will turn anEventSeries into a list of WhenEvents that NDSolve might handle better: eventQueue[ev_?TemporalDataEventSeriesQ, var_[t_]] := eventQueue[ev["Path"], var[t]]; eventQueue[bumps_List, var_[t_]] := WhenEvent[t == #1, var[t] -> var[t] + #2] & @@@ bumps;. (I do not know what limits there are on how many events NDSolve can handle efficiently.) – Michael E2 May 20 '15 at 14:25
• @Michael E2 - That is a nice idea, albeit it does get cumbersome. I probably will go for your Mod[t, 1] solution and adapt the EventSeries accordingly as stated. After all a model is a model and not 'reality' - some compromise has to be made eventually. – gwr May 20 '15 at 14:40
• Michael and Gwr, I have two questions, both based on the fact that the function models a real payment to an account : 1. Is there a way to substitute the t variable with a Date variable within the ODE? 2. If t=12 (or, depending on Q1, t=DateObject[{2020, 12}]`), how could I take into account within the ODE relation the double of the final year amount (13th month salary)? – Nate Mar 30 at 19:39