Programming a numerical method in the functional style

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 *)

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 $(x_o,y_o)$ 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 still easy to understand:

http://homepage.cem.itesm.mx/jose.luis.gomez/data/mathematica/metodos.pdf

Hope that helps