Classical Runge-Kutta, ODE4 | Solving ODEs in MATLAB

Here is the classical Runge-Kutta method. This was, by far and away, the world's most popular numerical method for over 100 years for hand computation in the first half of the 20th century, and then for computation on digital computers in the latter half of the 20th century. I suspect it's still in use today.

You evaluate the function four times per step, first in the beginning of the interval. And then use that to step into the middle of the interval, to get s2. Then you use s2 to step into the middle of the interval again.

And evaluate the function there again to get s3. And then use s3 to step clear across the interval, and get s4. And then take a combination of those four slopes, weighting the two in the middle more heavily, to take your final step. That's the classical Runge-Kutta method.

Here's our MATLAB implementation. And we will call it ODE4, because it evaluates to function four times per step. Same arguments, vector y out. Now we have four slopes-- s1 at the beginning, s2 halfway in the middle, s3 again in the middle, and then s4 at the right hand. 1/6 of s1, 1/3 of s2, 1/3 of s3, and 1/6 of s4 give you your final step. That's the classical Runge-Kutta method.

Carl Runge was a fairly prominent German mathematician and physicist, who published this method, along with several others, in 1895. He produced a number of other mathematical papers and was fairly well known.

Martin Kutta discovered this method independently and published it in 1901. He is not so nearly well known for anything else.

I'd like to pursue a simple model of combustion. Because the model has some important numerical properties. If you light a match, the ball of flame grows rapidly until it reaches a critical size. Then the remains at that size, because the amount of oxygen being consumed by the combustion in the interior of the ball balances the amount available through the surface.

Here's the dimensionless model. The match is a sphere, and y is its radius. The y cubed term is the volume of the sphere. And the y cubed accounts for the combustion in the interior.

The surface of the sphere is proportional y squared. And the y squared term accounts for the oxygen that's available through the surface. The critical parameter, the important parameter, is the initial radius, y0, y naught.

The radius starts at y0 and grows until the y cubed and y squared terms balance each other, at which point the rate of growth is 0. And the radius doesn't grow anymore. We integrate over a long time. We integrate over a time that's inversely proportional to the initial radius. That's the model.

Here's an animation. We're starting with a small flame here, a small spherical flame. You'll just see a small radius there. The time and the radius are shown at the top of the figure. It's beginning to grow.

When time gets to 50, we're halfway through. The flame sort of explodes, and then gets up the radius 1, at which time the two terms balance each other. And the flame stops growing. It's still growing slightly here, although you can't see it on this this scale.

Let's set this up for Runge-Kutta. The differential equation is y prime is y squared minus y cubed. Starting at zero, with the critical initial radius, I'm going to take to be 0.01. That means we're going to integrate out to two over y0 out to time 200.

I'm going to choose the step size to take 500 steps. I'm just going to pick that somewhat arbitrarily. OK, now I'm ready to use ODE4. And I'll store the results in y.

And it goes up to 1. I'm going to plot the results. So here's the values of t I need. And here's the plot.

Now, you can see the flame starts to grow. It grows rather slowly. And then halfway through the time interval, it's sort of explodes and goes up quickly, until it reaches a radius of 1, and then stays here.

Now this transition period is fairly narrow. And we're going to continue to study this problem. And is this transition area which is going to provide a challenge for the numerical methods.

Now here, we just went through it. We had a step size h, that we picked pretty arbitrarily. And we just generated these values. We have really little idea how accurate these numbers are.

They look OK. But how accurate are they? This is the critical question about the about the classical Runge-Kutta method. How reliable are the values we have here in our graph?

I have four exercises for your consideration. If the differential equation does not involve y, then this solution is just an integral. And the Runge-Kutta method becomes a classic method of numerical integration. If you've studied such methods, then you should be able to recognize this method.

Number. two-- find the exact solution of y prime equals 1 plus y squared, with y of 0 equals zero. And then see what happens with ODE4, when you try and solve it on the interval from t from 0 to 2.

Number three- what happens if the length of the interval is not exactly divisible by the step size? For example, if t final is pi, and the step size is 0.1. Don't try and fix this. It's just one of the hazards of a fixed step size.

And finally, exercise four-- investigate the flame problem with an initial radius of 1/1,000. For what value of t does the radius reach 90% of its final value?