Wolfram Cloud Document

CITE THIS NOTEBOOK: Imperative vs functional style code of Euler's method for ODEs by Eric Rimbey. Wofram Community DEC 28 2022.

Initialization and parameters

In[]:=

timeStep=.1;populationInit=10;timeInit=0;timeEnd=5;populationDerivative[pop_]:=.7pop

The following function is really the “gist” of the Euler method. If we know the derivative at a point, we can estimate the next point by extrapolating with the derivative. This was implicit in the Python code, but I want it to be explicit. Note that the output is another function, the extrapolation function for a given input value.

In[]:=

populationStepFn[pop_]:=Function[{dt},dtpopulationDerivative[pop]];

Imperative style

A hallmark of the imperative style is explicit management of state to execute algorithms, typically using iterator variables, counter variables, repeatedly updating result variables, etc. Mathematica has For, While, Do, and Until for looping in this imperative style. You might want to implement the algorithm with each of those constructs as a learning exercise. Here is a version using For.

In[]:=

For[(*initializationofstatevariables*)populationData={};population=populationInit;time=timeInit,(*terminalconditioncheck*)time<=timeEnd,(*incrementaftereachiteration*)population+=populationStepFn[population][timeStep];time+=timeStep,(*bodyofloop*)AppendTo[populationData,{time,population}]]

In[]:=

ForResult=populationData;(*savingtheresultforlatercomparison*)

Let’s look at the result

In[]:=

populationData//Short

Out[]//Short=

{{0,10},{0.1,10.7},48,{5.,294.57}}

Note that I didn’t create two separate lists, like the Python code, but instead created a single list that consisted of pairs. One reason is that I just like that better, these are points, so why separate them into separate lists for each coordinate? Another reason is that this is the form that the plotting function (ListPlot) will use for display.

Functional style

In the functional style, state management is tacitly handled by the implementation of the language. One typically defines functions to handle various logical elements of the algorithm, and ideally these functions have no side effects.

The approach here will be to start with an initial value and repeatedly apply a function. Mathematica uses the term “Nest” for this idea, and there are various Nest* functions. I’ll use NestWhileList. The basic construction looks like this in abstract:

In[]:=

NestWhileList[functionToBeRepeatedlyApplied,initialData,functionToDetermineWhenToStop]

There are additional arguments you can apply to guarantee termination and to add flexibility to how you determined done-ness.

The tricky bit is coming up with a “schema” for the data and a function that works on that data. I already know I want a list of points, so my initial data will be a 2-element list, and that should be based on our initial values defined above.

In[]:=

NestWhileList[functionToBeRepeatedlyApplied,{timeInit,populationInit},functionToDetermineWhenToStop]

Knowing when to stop is probably the next easiest thing. But remember, we’re not defining our own condition based on state, but a function that operates on the data as we iterate. As a condition, we want

currentTime<timeEnd

, so how to make that a function. Well, currentTime will be the first element of our most recently calculated point, and First is the function that extracts the first element of a list. There are built in functions for functionalizing the idea of comparison, and the one I want is LessThan. Composing these two gives us a nice, pure function we can apply to a single data point to do our iteration check.

In[]:=

NestWhileList[functionToBeRepeatedlyApplied,{timeInit,populationInit},LessThan[timeEnd]@*First]

Now we need a function that generates the next point from a previous point. So, for some input point

{time,population}

, we want to get

{time+timeStep,population+populationStepFn[population][timeStep]}

. So, we could do something like this:

In[]:=

NextData[{time_,population_}]:={time+timeStep,population+populationStepFn[population][timeStep]}

However, I’m going to use a different definition, on that’s still more functional. Notice that our above definition relied on populationStepFn and timeStep. But what if you want to improve the accuracy of Euler’s method by reducing the timeStep? You’d have to go back to the initialization section and change it’s value. To do any comparisons, this would get tedious. Similarly, what if you want to explore a slightly different function for the population derivative? Anyway, I’m going to add these as parameters to our NextData function, and I’m going to do it in a way that avoids us having to use anonymous function notation.

In[]:=

ClearAll[NextData];NextData[stepFn_,interval_][{time_,population_}]:={time+interval,population+stepFn[population][interval]}

Now we’ve got everything we need:

In[]:=

NestResult=NestWhileList[NextData[populationStepFn,timeStep],{timeInit,populationInit},LessThan[timeEnd]@*First];

In[]:=

NestResult//Short

Out[]//Short=

{{0,10},{0.1,10.7},48,{5.,294.57}}

Plotting

This is the easiest bit--we’ll just use ListPlot directly on our data.

In[]:=

ListPlot[ForResult]

Out[]=

In[]:=

ListPlot[NestResult]

Out[]=

Playing around

One nice advantage of the functional style is how easy it is to repeat things with just small tweaks. One reason for this is that functions evaluate to values (unlike looping procedures that just update state variables). Let’s say we want to refine our estimate, we can just “plug it in”:

In[]:=

ListPlot[NestWhileList[NextData[populationStepFn,.01],{timeInit,populationInit},LessThan[timeEnd]@*First]]

Out[]=

We can even extrapolate backward.

In[]:=

ListPlot[NestWhileList[NextData[populationStepFn,-.2],{timeInit,populationInit},GreaterThan[-3]@*First]]

Out[]=