Download solution of differntial equation using runge kutta method and more Lecture notes Mathematical Methods for Numerical Analysis and Optimization in PDF only on Docsity!
Euler’s Method, Taylor Series Method, Runge
Kutta Methods, Multi-Step Methods and Stability.
REVIEW: We start with the differential equation
dy(t) dt
= f (t, y(t)) (1.1) y(0) = y 0
This equation can be nonlinear, or even a system of nonlinear equations (in which case y is a vector and f is a vector of n different functions).
Numerical Solution of an ODE: The idea behind numerical solutions of a Differential Equation is to replace differentiation by differencing. A computer cannot differentiate but it can easily do a difference. (Differentiation is a continuous process. Differencing is a discrete process.) Now we introduce the most important tool that will be used in this section. By the time you’ve mastered this section, you’ll be able to do Taylor Expansions in your sleep. (I am already doing Taylor expansions in your sleep, right?!)
Taylor Series Expansion: You’ll recall (?) from your calculus class that if a function y(t) behaves nicely enough, then its Taylor series expansion converges:
y(t + ∆t) = y(t) + ∆ty′(t) +
∆t^2 y′′(t) +
∆t^3 y′′′(t) + ...
The Taylor series with remainder term is
y(t + ∆t) = y(t) + ∆ty′(t) +
∆t^2 y′′(t) +
∆t^3 y′′′(t) + ... +
n!
∆tny(n)(τ )
where τ is some value between t and t + ∆t. You can truncate this for any value of n.
Euler’s Method: If we truncate the Taylor series at the first term
y(t + ∆t) = y(t) + ∆ty′(t) +
∆t^2 y′′(τ ),
we can rearrange this and solve for y′(t)
y′(t) =
y(t + ∆t) − y(t) ∆t
Now we can attempt to solve (1.1) by replacing the derivative with a difference:
y((n + 1)∆t) ≈ y(n∆t) + ∆tf(n∆t, y(n∆t))
Start with y(0) and step forward to solve for any time.
What’s good about this? If the O term is something nice looking, this quantity decays with ∆t, so if we take ∆t smaller and smaller, this gets closer and closer to the real value. What can go wrong? The O term may be ugly. The errors can accumulate as I step forward
in time. Also, even though this may be a good approximation for y′(t) it may not converge to the right solution. To answer these questions, we look at this scheme in depth.
Terminology: From now on, we’ll call yn the numerical approximation to the solution y(n∆t); tn = n∆t. Euler’s method can then be written
yn+1 = yn + ∆tf(tn, yn) n = 1, ..., N − 1 (1.2)
This method assumes that you can move from one location to the next using the slope given by the equation (1.1). We saw last time that when we do this, our errors will decay linearly with ∆t. We will show this again today, but in two steps, so that we can generalize it. The proof should look very familiar!
Local Truncation Error: To be able to evaluate what we expect the order of a method to look like, we look at the
LT E(t) =
y(t + ∆t) − y(t) ∆t
− f(t, y(t)),
i.e. it is the residue when the exact solution of the ODE (1.1) is plugged into the numerical scheme. If yn is close to y(tn) then the LTE will be close to zero.
The local truncation error represents the terms neglected by truncating the Taylor series. This is not the error that we get from the method, (i.e. the difference between the real solution and the numerical solution) but will be connected.
If I don’t know y(t), what is the use of this definition? (and if I do know y(t), what do I need the method for?!). It turns out that even without explicit knowledge of the solution we can still calculate the LTE and use it as an estimate and control of the error, by placing certain smoothness assumptions on y(t) and using the Taylor Expansions.
Clearly, at time tn, Euler’s method has Local Truncation Error:
LT E =
y(tn + ∆t) − y(tn) ∆t
− f(tn, y(tn)) = O(∆t),
in other words, we can write this
y(tn+1) = y(tn) + ∆tf(tn, y(tn)) + ∆tLT E.
Of course, the method is yn+1 = y(tn) + ∆tf(tn, yn).
Subtract these two,
|y(tn+1) − yn+1| = |y(tn) − yn + ∆t (f(tn, y(tn)) − f(tn, yn)) + ∆tLT E| ≤ |y(tn) − yn| + ∆t |f(tn, y(tn)) − f(tn, yn)| + ∆t|LT E| ≤ |y(tn) − yn| + ∆tL |y(tn) − yn| + ∆t|LT E|.
Because f is Lipschitz continuous,
f(tn, y(tn)) − f(tn, yn) y(tn) − yn
| ≤ L.
Taylor Series Methods: To derive these methods we start with a Taylor Expansion:
y(t + ∆t) ≈ y(t) + ∆ty′(t) +
∆t^2 y′′(t) + ... +
r!
y(r)(t)∆tr.
Let’s say we want to truncate this at the second derivative and base a method on that. The scheme is, then:
yn+1 = yn + fn∆t +
f t′n 2
∆t^2.
The Taylor series method can be written as
yn+1 = yn + ∆tF (tn, yn, ∆t)
where F = f + 12 ∆tf′. If we take the LTE for this scheme, we get (as expected)
LT E(t) =
y(tn + ∆t) − y(tn) ∆t
− f(tn, y(tn)) −
∆tf′(tn, y(tn)) = O(∆t^2 ).
Of course, we designed this method to give us this order, so it shouldn’t be a surprise!
So the LTE is reasonable, but what about the global error? Just as in the Euler Forward case, we can show that the global error is of the same order as the LTE. How do we do this? We have two facts, y(tn+1) = y(tn) + ∆tF (tn, y(tn), ∆t),
and yn+1 = yn + ∆tF (tn, yn, ∆t)
where F = f + 12 ∆tf′. Now we subtract these two
|y(tn+1) − yn+1| = |y(tn) − yn + ∆t (F (tn, y(tn)) − F (tn, yn)) + ∆tLT E| ≤ |y(tn) − yn| + ∆t |F (tn, y(tn)) − F (tn, yn)| + ∆t|LT E|.
Now, if F is Lipschitz continuous, we can say
en+1 ≤ (1 + ∆tL)en + ∆t|LT E|.
Of course, this is the same proof as for Euler’s method, except that now we are looking at F , not f, and the LT E is of higher order. We can do this no matter which Taylor series method we use, how many terms we go forward before we truncate.
Advantages and Disadvantages of the Taylor Series Method:
advantages a) One step, explicit b) can be high order c) easy to show that global error is the same order as LTE disadvantages Needs the explicit form of derivatives of f.
Runge-Kutta Methods To avoid the disadvantage of the Taylor series method, we can use Runge-Kutta methods. These are still one step methods, but they depend on estimates of the solution at different points. They are written out so that they don’t look messy:
Second Order Runge-Kutta Methods:
k 1 = ∆tf(ti, yi) k 2 = ∆tf(ti + α∆t, yi + βk 1 ) yi+1 = yi + ak 1 + bk 2
let’s see how we can chose the parameters a,b, α, β so that this method has the highest order LT E possible. Take the Taylor expansions to express the LTE:
k 1 (t) = ∆tf(t, y(t)) k 2 (t) = ∆tf(t + α∆t, y + βk 1 (t) = ∆t
( f(t, y(t) + ft(t, y(t))α∆t + fy (t, y(t))βk 1 (t) + O(∆t^2 )
)
LT E(t) =
y(t + ∆t) − y(t) ∆t
a ∆t
f(t, y(t))∆t −
b ∆t
(ft(t, y(t))α∆t + fy (t, y(t)βk 1 (t)
=
y(t + ∆t) − y(t) ∆t
− af(t, y(t)) − bf(t, y(t)) − bft(t, y(t))α
− bfy(t, y(t)βf(t, y(t)) + O(∆t^2 )
= y′(t) +
∆ty′′(t) − (a + b)f(t, y(t)) − ∆t(bαft(t, y(t)) + bβf(t, y(t))fy(t, y(t)) + O(∆t^2 )
= (1 − a − b)f + (
− bα)∆tft + (
− bβ)∆tfyf + O(∆t^2 )
So we want a = 1 − b, α = β = (^21) b.
Fourth Order Runge-Kutta Methods:
k 1 = ∆tf(ti, yi) (1.3)
k 2 = ∆tf(ti +
∆t, yi +
k 1 ) (1.4)
k 3 = ∆tf(ti +
∆t, yi +
k 2 ) (1.5) k 4 = ∆tf(ti + ∆t, yi + k 3 ) (1.6)
yi+1 = yi +
(k 1 + k 2 + k 3 + k 4 ) (1.7)
The second order method requires 2 evaluations of f at every timestep, the fourth order method requires 4 evaluations of f at every timestep. In general: For an rth order Runge- Kutta method we need S(r) evaluations of f for each timestep, where
S(r) =
r for r ≤ 4 r + 1 for r = 5 and r = 6 ≥ r + 2 for r ≥ 7
