Numerical Methods - Runge-Kutta Method
Numerical Methods for Solving Differential Equations
The Runge-Kutta Method
Theoretical Introduction
In the last lab you learned to use Heuns's Method to generate a numerical solution to an initial value problem of the form:
²â′Ìý=Ìýf(³æ,Ìý²â)
y(xo)Ìý=Ìýyo
We discussed the fact that Heun's Method was an improvement over the rather simple Euler Method, and that though it uses Euler's method as a basis, it goes beyond it, attempting to compensate for the Euler Method's failure to take the curvature of the solution curve into account. Heun's Method is one of the simplest of a class of methods calledÌýpredictor-corrector algorithms. In this lab we will address one of the most powerful predictor-corrector algorithms of all—one which is so accurate, that most computer packages designed to find numerical solutions for differential equations will use it by default—theÌýfourth order Runge-Kutta Method. (For simplicity of language we will refer to the method as simply theÌýRunge-Kutta MethodÌýin this lab, but you should be aware that Runge-Kutta methods are actually a general class of algorithms, the fourth order method being the most popular.)
The Runge-Kutta algorithm may be very crudely described as "Heun's Method on steroids." It takes to extremes the idea of correcting the predicted value of the next solution point in the numerical solution.Ìý(It should be noted here that the actual, formal derivation of the Runge-Kutta Method willÌýnotÌýbe covered in this course. The calculations involved are complicated, and rightly belong in a more advanced course in differential equations, or numerical methods.)ÌýWithout further ado, using the same notation as in the previous two labs, here is a summary of the method:
Even though we do not plan on deriving these formulas formally, it is valuable to understand the geometric reasoning that supports them. Let's briefly discuss the components in the algorithm above.
First we note that, just as with the previous two methods, the Runge-Kutta method iterates theÌýx-values by simply adding a fixed step-size ofÌýhÌýat each iteration.
TheÌýy-iteration formula is far more interesting. It is a weighted average of four values—k1,Ìýk2,Ìýk3, andÌýk4. Visualize distributing the factor of 1/6 from the front of the sum. Doing this we see thatÌýk1ÌýandÌýk4Ìýare given a weight of 1/6 in the weighted average, whereasÌýk2ÌýandÌýk3Ìýare weighted 1/3, or twice as heavily asÌýk1ÌýandÌýk4. (As usual with a weighted average, the sum of the weights 1/6, 1/3, 1/3 and 1/6 is 1.) So what are theseÌýkiÌývalues that are being used in the weighted average?
k1Ìýyou may recognize, as we've used this same quantity on both of the previous labs. This quantity,ÌýhÌýf(xn,Ìýyn), is simply Euler's prediction for what we've previously called Δy—the vertical jump from the current point to the next Euler-predicted point along the numerical solution.
k2Ìýwe have never seen before. Notice theÌýx-value at which it is evaluating the functionÌýf.ÌýxnÌý+Ìýh/2 lies halfway across the prediction interval. What about theÌýy-value that is coupled with it?ÌýynÌý+Ìýk1/2 is the currentÌýy-value plus half of the Euler-predicted ΔyÌýthat we just discussed as being the meaning ofÌýk1. So this too is a halfway value, this time vertically halfway up from the current point to the Euler-predicted next point. To summarize, then, the functionÌýfÌýis being evaluated at a point that lies halfway between the current point and the Euler-predicted next point. Recalling that the functionÌýfÌýgives us theÌýslopeÌýof the solution curve, we can see that evaluating it at the halfway point just described, i.e.Ìýf(xnÌý+Ìýh/2,ÌýynÌý+Ìýk1/2), gives us an estimate of the slope of the solution curve at this halfway point. Multiplying this slope byÌýh, just as with the Euler Method before, produces a prediction of theÌýy-jump made by the actual solution across the whole width of the interval, only this time the predicted jump is not based on the slope of the solution at the left end of the interval, but on the estimated slope halfway to the Euler-predicted next point. Whew! Maybe that could use a second reading for it to sink in!
k3Ìýhas a formula which is quite similar to that ofÌýk2, except that whereÌýk1Ìýused to be, there is now aÌýk2. Essentially, theÌýf-value here is yet another estimate of the slope of the solution at the "midpoint" of the prediction interval. This time, however, theÌýy-value of the midpoint is not based on Euler's prediction, but on theÌýy-jump predicted already withÌýk2. Once again, this slope-estimate is multiplied byÌýh, giving us yet another estimate of theÌýy-jump made by the actual solution across the whole width of the interval.
k4ÌýevaluatesÌýfÌýatÌýxnÌý+Ìýh, which is at the extreme right of the prediction interval. TheÌýy-value coupled with this,ÌýynÌý+Ìýk3, is an estimate of theÌýy-value at the right end of the interval, based on theÌýy-jump just predicted byÌýk3. TheÌýf-value thus found is once again multiplied byÌýh, just as with the three previousÌýki, giving us a final estimate of theÌýy-jump made by the actual solution across the whole width of the interval.
In summary, then, each of theÌýkiÌýgives us an estimate of the size of theÌýy-jump made by the actual solution across the whole width of the interval. The first one uses Euler's Method, the next two use estimates of the slope of the solution at the midpoint, and the last one uses an estimate of the slope at the right end-point. EachÌýkiÌýuses the earlierÌýkiÌýas a basis for its prediction of theÌýy-jump.
Ìý
This means that the Runge-Kutta formula forÌýyn+1, namely:
yn+1Ìý=ÌýynÌý+Ìý(1/6)(k1Ìý+Ìý2k2Ìý+Ìý2k3Ìý+Ìýk4)
is simply theÌýy-value of the current point plus a weighted average of four differentÌýy-jump estimates for the interval, with the estimates based on the slope at the midpoint being weighted twice as heavily as the those using the slope at the end-points.
As we have just seen, the Runge-Kutta algorithm is a little hard to follow even when one only considers it from a geometric point of view. In reality the formula was not originally derived in this fashion, but with a purely analytical approach. After all, among other things, our geometric "explanation" doesn't even account for the weights that were used. If you're feeling ambitious, a little research through a decent mathematics library should yield a detailed analysis of the derivation of the method.
It's now time toÌýimplement the method inÌýMathematica.
 |
If you're lost, impatient, want an overview of this laboratory assignment, or maybe even all three, you can click on the compass button on the left to go to the table of contents for this laboratory assignment. |