So previously we've looked at some simple ways to solve differential equations numerically,
specifically first-order ordinary differential equations. We're going to get a little deeper
now and look at the standard bread and butter way to solve these equations. It's a little
more complicated, but it's a lot more accurate, it other ways we need relatively few function
evaluations in order to calculate the value of the function's differential equation we're
given. This particular method is the Runge-Kutta 4th order method, there's actually a family
of Runge-Kutta 4th order methods, they all come from writing out the Taylor series expansion with
the Runge-Kutta formula, and solving for the constraints. There's a number of possible
formulas it could write out. This is the one that people generally use when they talk about
a Runge-Kutta 4th order algorithm because it's particularly simple to write out, the
constants are relatively easy to remember. So we're going to start out with our definition
of our differential equation that we're trying to solve. We have dx/dt is some function of
x and t, and we have an initial point at which we know the value of our function. What we
want to do is solve for x as a function of t for all t, and what we're going to do is
find a way to figure out what x is at some value of t(0)+dt, because if we can do that
we can keep repeating that on and on until we have a function x of t at all time t. So
all the Runge-Kutta formulas are of this form. x(t)+dt is equal to x of the current value
of t plus some interval dt, and the average slope over that interval, and that's where
the difference comes in, different Runge-Kutta methods have different ways of calculating
the slope. This particularly version that we're going to look at the average slope is
calculated as one sixth of four slopes, it's 4th order, and to get it into 4th order you
need to evaluate the function, i.e. the function f(x,t), which is dx/dt, in four places, and
then you weight that average, so in this case we take the slope of the first point twice
the slope of the second point, twice the slope of the third point, plus the slope at the
fourth point, and then we divide by six to average it out. You can tell that it's six
because you have 1+2+2+1 total slopes to look at, that adds to sixth, so if we want to have
something that's the same order of magnitude as the slope we're going to divide by the
sum of the weights which is six, we divide by six, and that's our average slope, so that
part's pretty easy to memorize, you just remember the two in the middle are doubled, and then
we have to divide by six. Now the next question is at what points x and t do we evaluate these
four different derivatives. Well, we start out at our initial time t(0), and that's the
time at which we know what the value of x is, so we have that point. We want to find x as
a function of all t, we start out with only one t. Clearly this is the only point we know,
this is the point we have to start out with when we calculate our derivative, so that's
our k(1), it's the function evaluated at x(0), and t(0), in other words dx/dt evaluated at
that point, so if we know the slope of the function exactly at that point it's an approximation
to the derivative of the function over the next time period. What we're going to do to
get the next point is go forward a half step, dt/2, assuming that the function is roughly
linear over that time step, and calculate an estimate for the function at t(0)+dt/2,
it's not a perfect value of the function because we're assuming it's linear, but it's a good
estimate. That's our first estimate for the function ,and so now we want to calculate
the derivative of the function at that point, well that's simple enough, if we just take
x, which is x(0)+dt/2 *k(1) for our x, and then for our t value its t(0)+dt/2. So we
calculate our derivative of the function evaluated at that point. So that's our second slope.
This should be a better estimate of the slope, the midpoint, than the value of the slope
at the first point of the interval, so now we go from that point over and we take a half
step using our second slope. So now we have a third point, and that third point is just
going to be x(0)+k(2)*dt/2, that's the value of x, and the value of t is of course t(0)+dt/2,
and we can calculate our slope there, or in other words the derivative of the function
evaluated at that point, and that gives us a new estimate for the derivative at the midpoint,
so then we take that estimate of the slope, the estimate of the derivative of the function,
go back to our first point, and take a full step, so we take a full step over for the
full dt interval and get a point there. So that gives us our fourth estimated point
for the function, and we can calculate our slope at that fourth point, the fourth point
which is x(0)+K(3)*dt, a full dt now, not just half a step, and then the t
value is t(0)+dt, so that gives us our fourth slope, our fourth value of the derivative
of the function at four different estimated points. Now for a real function these would
all be much closer together, if we want to have a good estimate of the slope we have
to have a very small dt or else the function would change too much, so for a real
function this would clearly be a dt that's too big, because the derivative is
changing drastically over the course of the time step, but we're exaggerating everything,
so you can see the differences in the slopes at the different points better. That gives
us the information we need. OK, so now we can go through and we can see the average
slope is one-sixth of one times the first slope plus two time the second slope, plus
two times the third slope, plus the fourth slope, we add those together. So what does
that mean, what does it mean when we add together a line? Well really each line is just the
value of the slope, it's just a number, and so numerically, well, the first one looks about
1.2, the second one maybe 0.5, the third one is flatter, so maybe 0.3, and the last one
is, it's higher than the first one so we'll call it 2, and we're just estimating here,
so we average those together and that would give us a slope of about 0.8. We computed
the average slope, now we can make the full step and give us our Runge-Kutta estimate
of what the value of the function is at t(0)+dt. We're going to just minimize that little window over
in the corner so we can be remembering that as we go through and plug in the numbers for
this, so again remember this is what we're doing, we're computing an average slope,
multiplying that by dt, and then incrementing our current value with that average slope,
and let's use the same function we used for our lower order solution algorithms, it's
3*x^(0.5)+t^3, so this is a function of both x and t. And of course we need to have a time
step size, so we'll choose 0.1, and we need to have an initial point, so we're going to
start at t(0.2), and we're going to assume that it's 0.1. Now, these number do not necessarily
correspond directly to the numbers on the graph, don't try to look at the graph and
get the numbers to fit exactly. In order to fit everything on the graph I had to choose
a bit of a different scaling for that, so just use the figure as a reference to remember
what goes where, not to actually look up what the value of the slope is of the function
at any given point, that's just to give you a sense of what order the function values
are. Now our first slope in red we evaluate from the first point, so we go through, we
plug in x(0.2), that's our initial point, and t=0.2. Well x(0.2) is just 0.1, so we
plug that in for x, and we get a slope of 0.95668, great, that's our first slope, so
now we're going to calculate the value of the function at t+dtd/2.dt is
0.1, so dt/2 is 0.05, so now we're interested in the value of the function at 0.25, and
so that's x(0.2)+k(1)*dt/2, and we get our first estimate of x(1) at t=0.25, so now
we can figure out the value of the slope, that's the light blue line, at that point,
we just plug in the numbers, we can plug it in using the formula we're given, we get an
estimate of slope is 1.16912, that's not that far from 0.95688, so we probably have
a decent time step, this slope isn't changing a huge amount. So now we can figure out our
2nd guess for the value of the function at the midpoint, and so we take the initial value
plus k(2)*dt/2, plug that in we get a new value of the function, and from that
we can figure out the slope at that new estimate of the function x at the midpoint 0.25, plug
it in we get a new slope, it's 1.20982, the slopes aren't changing very much, so that's
good, we've got a decent estimate, and then finally we take a full time step, we have
the value at 0.02+k(3), the orange/yellow slope times dt we make our full time
step, and we get a value of 0.22098, OK, great, and we can calculate the slope at that point,
plug in all our numbers, and we get a slope of 1.43726. Sounds reasonable, slope hasn't
changed too much. So we got our four slopes, and now we just average them together, and
we get an average slope over the entire interval dt of 1.19170, so as we see that's,
you know, intermediate between all these different values, so that seems like it could be a pretty
good estimate, it doesn't look like we're taking too big a times step. So our final
answer, now we've got our average slope, it's pretty easy. We just say it's going to be
x(0), which we know, plus the average slope times dt, plug that in, this math is
easy, we're just plugging that in, and we get a value of 0.21917. Is that good? Well,
what we can do is actually go through and use ODE of four or five, set a really tight
tolerance, in this case I've set the tolerance down to 10^-9, so it's a very accurate answer,
we get the true answer is 0.219257. So we got very close to that true answer using just
one time step and four function evaluations. So if we look at that we're within 0.04% of
the true answer with just one times step over that period, so it's a very accurate method
for most differential equations. You could come up with stiff differential equations,
or some other sort of quirky cases, which do come up from time to time, where you want
to try something else, but this ends up working very efficiently in most cases So that is
our summary of a particular example, the most common example, the one that's sort of the
default that you try, of 4th order Runge-Kutta algorithms to calculate numerical solutions
to first-order ordinary differential equations.