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["DifferentialEquations`NDSolveUtilities`"];
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.
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 letNDSolve
choose its own method, unless after doing so, it seems to give bogus results. $\endgroup$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. $\endgroup$b
andg
as follows:b[t_?NumericQ]:=pulseTrain[0,1/32,1,10][t];g[t_?NumericQ]:=32 b[t]
$\endgroup$"DiscontinuityProcessing" -> False
with a step size of 1/32, which is what happens when you have function protected byNumericQ
. TheDiracDelta
was a gratuitous addendum, which cannot work with regular Euler scheme, without approximating it by a finite pulse, as you do with yourb[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. $\endgroup$