Nail-Tinted Glasses

3rd order Runge-Kutta

Now, when I come across something like a Runge Kutta formulation, it tends to go into my notebook. Once upon a time, I was capable of doing math, and I'm trying to reacquire the skill. I don't yet feel the need to work it out for myself, as much as go through the derivation exercise once, to force myself to understand it.

In other words, I hit the web, looking for an explanation.

Searching reveals a number of explanations of the general approach, several derivations of 2nd order RK (because that's the easy one to derive), and discussion of 4th order RK (because that's the useful one, and the formula is pretty, if you are into that sort of thing). But almost nothing on 3rd order, and nothing to help somebody like me who gets the idea, but can't get the symbols to come out right.

Once I finish this, there will be.

Runge-Kutta is a recipe for finding a numerical solution to an ordinary differential equation. An ordinary differential equation is one with a single independent variable ( y' = F(x,y) ), and numerical analysis is frequently used for problems that don't offer an exact solution. For example, you might use one for a many body problem, trying to track where the planets are as a function of time.

With nice functions (continuous, differentiable, low sodium), and small steps in time, the next value is going to be relatively close to the current value, and if you cleverly choose a simple function that is similar to the real one, you can get a good approximation of the next value by using the next value of the simple function.

A Taylor series approximation is great at this sort of thing. Depending on how many derivatives you have available, you can approximate the actual function with a polynomial, and work from there.
y(x+h) = y(x) + hy'(x) + h^2y''(x)/2! + h^3y'''(x)/3! + ...

Unfortunately, in this case we have only one derivative (y'=F) to work with. If we want an answer that's better than the first order approximation
y(x+h) = y(x) + hy'(x)
then we are going to have to work for it.

At a high level, what RK does is find a weighted average for the slope. That is, instead of accepting the slope at the origin as gospel, the method samples a few other slopes that might be along the way, and uses a weighted average. The meat of the derivation is in the weightings - the goal is to choose weightings that agree with the Taylor series expansion of the appropriate order of h.

The general form of an RK solution is

Y[n+1] = Y[n] + h( w[1]k[1] + w[2]k[2] + w[3]k[3] + ... )

where w[j] is just a weight, and k[j] is a function of F and the previous k[] expressions.

k[j] = F( x[n] +ha[j], Y[n]+ h Sum B[ji]k[i] )

So we'll need to be able to evaluate F(x+a, y+b) - we can use a Taylor expansion for that. For second order RK derivations you only need first order terms of the Taylor series. But because we are here deriving a third order formulation, we need to keep all terms of h^3. There's a factor h out in front of all of the wk's, so we can discard the high order terms, but we do need to work the expansion out to second order.

F(x+a,y+b) = F + aFx + bFy + 1/2 ( a^2Fxx + 2abFxFy + b^2Fyy )

where Fx and Fy are the partial derivatives in those directions, and everything is evalated at (x,y). In practice, we aren't going to solve those partial derivatives - we're just collecting them with weights in front of them, so that at the end of the exercise we can choose the weights to match the Taylor series of our target function.

It seems that convention dictates that a[1] = 0, and therefore

k[1] = F
 
k[2] 
= f( x + ha[2], Y + hB[21]k[1] )
= F + h (a[2]Fx + B[21]FFy) + h^2/2 ( a[2]^2Fxx + 2a[2]B[21]FFxy + B[21]^2F^2Fyy )
 
k[3]
= f( x + ha[3], Y + h(B[31]k[1] + B[32]k[2]) )
 
= F + h (a[3]Fx + (B[31]F + B[32]k[2])Fy) 
+ h^2/2 ( a[3]^2Fxx + 2a[3](B[31]F+B[32]k[2])Fxy + (B[31]F+B[32]k[2])^2Fyy)
 
= F + h (a[3]Fx + B[31]FFy)
+ h B[32]k[2]Fy
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h^2/2 ( 2a[3]B[32]k[2]Fxy + 2B[31]FB[32]k[2]Fyy +  B[32]^2k[2]^2Fyy )
 
= F + h (a[3]Fx + B[31]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h B[32]Fyk[2]
+ h^2/2 ( 2a[3]B[32]Fxy + 2B[31]FB[32]Fyy ) k[2]
+ h^2/2 B[32]^2Fyy k[2]^2

We are looking for a solution to third order in h. There's a factor of h waiting for us at the end, so we only need second order here. Thus, we can take a low order approximation of k[2].

= F + h (a[3]Fx + B[31]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h B[32]Fy ( F + h(a[2]Fx + B[21]FFy) + O(h^2) )
+ h^2/2 ( 2a[3]B[32]Fxy + 2B[31]FB[32]Fyy ) ( F + O(h) )
+ h^2/2 B[32]^2Fyy (F + O(h))^2
 
= F + h (a[3]Fx + B[31]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h B[32]Fy ( F + h(a[2]Fx + B[21]FFy ))
+ h^2/2 ( 2a[3]B[32]Fxy + 2B[31]FB[32]Fyy )F
+ h^2/2 B[32]^2FyyF^2
 
= F + h (a[3]Fx + B[31]FFy + B[32]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h^2 B[32]Fy ( a[2]Fx + B[21]FFy )
+ h^2/2 ( 2a[3]B[32]Fxy + 2B[31]FB[32]Fyy )F
+ h^2/2 B[32]^2FyyF^2
 
= F + h (a[3]Fx + B[31]FFy + B[32]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy + B[31]^2F^2Fyy )
+ h^2 B[32]Fy ( a[2]Fx + B[21]FFy )
+ h^2/2 ( 2a[3]B[32]FFxy + 2B[31]B[32]F^2Fyy )
+ h^2/2 B[32]^2FyyF^2
 
= F + h (a[3]Fx + B[31]FFy + B[32]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy  )
+ h^2 B[32]Fy ( a[2]Fx + B[21]FFy )
+ h^2/2 ( 2a[3]B[32]FFxy  )
+ h^2/2 ( B[32]^2FyyF^2 + 2B[31]B[32]F^2Fyy + B[31]^2F^2Fyy )
 
= F + h (a[3]Fx + B[31]FFy + B[32]FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3]B[31]FFxy  )
+ h^2 B[32]Fy ( a[2]Fx + B[21]FFy )
+ h^2/2 ( 2a[3]B[32]FFxy  )
+ h^2/2 ( B[31] + B[32] ) ^ 2 F^2 Fyy
 
= F + h (a[3]Fx + ( B[31] + B[32] )FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3](B[31]+B[32])FFxy  )
+ h^2 B[32]Fy ( a[2]Fx + B[21]FFy )
+ h^2/2 ( B[31] + B[32] ) ^ 2 F^2 Fyy
 
= F + h (a[3]Fx + ( B[31] + B[32] )FFy)
+ h^2/2 ( a[3]^2Fxx + 2a[3](B[31]+B[32])FFxy  )
+ h^2 ( a[2]B[32]FxFy  + B[32]B[21]FFy^2 )
+ h^2/2 ( B[31] + B[32] ) ^ 2 F^2 Fyy
 
 
Y[n+1] 
= Y 
+ hF               (w1+w2+w3)
+ h^2Fx            (w3 a[3]              + w2 a[2])
+ h^2FFy           (w3 (B[31]+B[32])     + w2 B[21] )
+ h^3/2 * Fxx      (w3 a[3]^2            + w2 a[2]^2 )
+ h^3/2 * FFxy   2 (w3 a[3](B[31]+B[32]) + w2 a[2]B[21] )
+ h^3/2 * FFFyy    (w3 (B[31]+B[32])^2   + w2 B[21]^2 )
+ h^3/2 * FxFy   2 (w3 a[2]B[32])
+ h^3/2 * FFyFy  2 (w3 B[32]B[21])
 
Y[n+1] 
= Y + hY' + (h^2/2)Y'' + (h^3/6)Y'''
= Y + hF + (h^2/2)F' + (h^3/6)F''
= Y + hF + (h^2/2)(Fx+FyF) + (h^3/6)(d/dx)(Fx+FyF)
= Y + hF + (h^2/2)(Fx+FyF) + (h^3/6)(Fxx + FxyF + FxyF + FyyFF + FyFx + FyFyF)
= Y + hF + (h^2/2)(Fx+FyF) + (h^3/6)(Fxx + 2 FFxy  + FFFyy + FxFy + FFyFy)

Now we just choose coefficients so that each term of these two expression is equal.

w1 = 1/6
w2 = 4/6
w3 = 1/6
a[2] = 1/2
a[3] = 1
B[21] = 1/2
B[31] = -1
B[32] = 2
 
k[1] = f(x,Y)
k[2] = f(x+h/2,Y+k[1]/2)
k[3] = f(x+h, Y-k[1]+2k[2])

January 18, 2004 11:03 PM | TrackBack