# Programming a numerical method in the functional style

I am new to Mathematica and I would like to learn a bit more about functional programming.

At the moment I have assignments like programming different numerical methods (for integration: trapezoidal rule, simpson's rule ,.. , for differential equations: euler's method, midpoint Runge-Kutta, ..). I've been implementing all of the methods with procedural programing, with While loops and For loops, but now with Eluer's method I've started using functional programming and have used Table and FoldList instead of a For loop. It comes together quite nicely. Here is an example:

euler[f_Symbol, y0_, {tz_, tk_}, h0_: 0.1] :=
Module[{t0 = tz, h = h0, n = Abs[tz - tk]/h0, Y, T, y},
y = FoldList[#1 + h*f[#1, #2] &, y0, Table[t0 + i*h, {i, 0, n - 1}]];
T = Table[t0 + i*h, {i, 0, n + 1}];
Y = Table[N[{T[[j]], y[[j]]}], {j, 1, n + 1}];
ListPlot[Y, Joined -> True, PlotRange -> All]
]


I would like to do something similar with the midpoint Runge-Kutta method, but I don't know how. Here is the procedural code:

midpointRK[f_, y0_, ta_, tb_, h0_: 0.1] :=
Module[{tA = ta, tB = tb, h = h0, n, Y, g, y},
n = Abs[tB - tA]/h;
y[0] = y0;
For[i = 0, i <= n, i++,
y[i + 0.5] = y[i] + (h/2)*f[y[i], tA + i*h];
g = f[y[i + 0.5], tA + i*h*0.5];
y[i + 1] = y[i] + h*g;
Y = Table[{N[tA + j*h], y[j]}, {j, 0, n}]
];
ListPlot[Y, Joined -> True, PlotRange -> All]
]


Can someone please help me in this regard and explain to me a little, what each command does?

-

Since you're learning, I'll show you how to break down your procedural function midpointRK into bare parts and reassemble it in a functional style. FoldList will again be the function of choice here, so since you're already familiar with it, I'll skip the explanation on that.

First, observe that ta, tb and h0 are all constants, so there is no need to reassign them — you can use them as is. n, f and y0 are also constants, so if you look at the body of your For loop, it requires

• The current iteration step, i
• A current state (or value), y[i]

and computes

• An intermediate value y[i + 0.5] from y[i] and the other constants
• An intermediate value g from y[i + 0.5] and other known constants
• An "updated" state, y[i + 1] (for the next iteration) from y[i], g and other known constants

The Y is being overwritten each time and not used, so it makes sense to include it only at the end. So since your intermediate values are only being used in the current iteration and not in future iterations, and you only care about y[i] for integer i, it really could be stored in any variable (not necessarily y[i + 0.5]).

We can now encapsulate what we have so far in a little wrapper function that takes the known constants, the current state and current iteration index and outputs the updated state:

rk[ta_, h_, n_, f_][y_, i_] := With[{ymid = y + h f[y, ta + i h]/2},
y + h f[ymid, ta + i h/2]]


The above should be self-explanatory. Now we can package this up into a FoldList and "fold" over n iteration steps as:

midpointRKFunctional[f_, y0_, ta_, tb_, h0_: 0.1] :=
With[{n = Abs[tb - ta]/h0},
Transpose@{ta + h0 Range[0, n], FoldList[rk[ta, h0, n, f], y0, Range[0, n - 1]]}
]


The final output from this function is the same as the Y in your procedural code. You can verify that it indeed does return the same result:

midpointRK[#^2 &, 0.1, 2, 1, 0.1] == midpointRKFunctional[#^2 &, 0.1, 2, 1, 0.1]
(* True *)

-
Perhaps this can be useful if aren't familiar with With[] and companions ... –  belisarius Mar 31 '13 at 14:52
Thank you very much for your help, it worked nicely. –  Gregor Mar 31 '13 at 15:11
P.S. For those wondering, I modified the OP's midpointRK to return Y instead of a plot. –  The Toad Mar 31 '13 at 18:05

This is an implementation of Eulers method i put together using NestList.

NestList and Euler's method

Hope that helps. I would also recommend Paul Wellins book called Programming with Mathematica. I am new to mathematica but this book has helped me a lot. I am by no means an expert but to really harness the functional style you need to understand things like Nest, NestList,Map, Thread etc etc. I would start with simple example first and just get familiar with these features and then apply them to numerical methods.

I found this site useful as well, it doesn't use many of the functional constructs but its still useful to see how you would implement using more traditional structured programming methods.

http://math.fullerton.edu/mathews/numerical.html

David.

-

If you use NestWhileList together with Function in order to implement numerical methods, the first argument of NestWhileList is the numerical method, the second argument is the initial point, and the third argument is the condition to continue calculating. Furthermore, you get as an output the list of all the points that make your numerical solution. This is an implementation of Euler's method (the equation to solve is of the form y'=f(x,y), the initial point is (xo,yo) and h is the step size):

euler[f_, {xo_, yo_}, xend_, h_] :=
NestWhileList[
Function[pointxy, pointxy + {h, h*Apply[f, pointxy]}],
{xo, yo},
Function[pointxy, pointxy[[1]] < xend]]


After evaluating the implementation, you can use it like this:

f[x_, y_] := 0.1*y*(100 - y);
data= euler[f, {0, 50}, 0.5, 0.02];
ListPlot[data]


Some options for ListPlot improve the output:

ListPlot[data, AxesOrigin -> {0, 0}, PlotLabel -> y'[x] == f[x, y[x]] ]


Here is the link if you want to see a tutorial going from a basic use of NestWhileList to its use in Euler, Improved Euler, Higher order Euler, etc. The comments and titles are in my language (Spanish), however I believe the examples are anyway easy to understand:
http://homepage.cem.itesm.mx/jose.luis.gomez/data/mathematica/metodos.pdf

Hope that helps

-