Study with the several resources on Docsity
Earn points by helping other students or get them with a premium plan
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Material Type: Notes; Professor: Datta; Class: Numerical Analysis; Subject: MATHEMATICAL SCIENCES; University: Northern Illinois University; Term: Spring 2009;
Typology: Study notes
1 / 24
This page cannot be seen from the preview
Don't miss anything!
We consider the problem of numerically solving a system of differential equations of the form dy dt =^ f^ (t, y),^ a^ ≤^ t^ ≤^ b,^ y(a) =^ α^ (given)^. Such a problem is called the Initial Value Problem or in short IVP, because the initial value of the solution y(a) = α is given.
Since there are infinitely many values between a and b, we will only be concerned here to find approximations of the solution y(t) at several specified values of t in [a, b], rather than finding y(t) at every value between a and b.
Denote
t 0 = a < t 1 < t 2 < · · · tN = b. ^
PSfrag replacements
a = t 0 t 1 t 2 tN = b
Theorem: (Existence and Uniqueness Theorem for the IVP). The initial value problem: (^) { y′^ = f (t, y) y(a) = α has a unique solution y(t) for a ≤ t ≤ b, if f (t, y) is continuous on the domain, given by R = {a ≤ t ≤ b, ∞ < y < ∞} and satisfies the following inequality: |f (t, y) − f (t, y∗)| ≤ L|y − y∗|,
Whenever (t, y) and (t, y∗) ∈ R. Definition. The condition |f (t, y) − f (t, y∗)| ≤ L|y − y∗| is called the Lipschitz Condi- tion. The number L is called a Lipschitz Constant.
Definition. A set S is said to be convex if whenever (t 1 , y 1 ) and (t 2 , y 2 ) belong to S, the point ((1 − λ)t 1 + λt 2 , (1 − λ)y 1 + λy 2 ) also belongs to S for each λ when 0 ≤ λ ≤ 1.
Simplification of the Lipschitz Condition for the Convex Domain If the domain happens to be a convex set, then the condition of the above Theorem reduces to (^) ∣ ∣∣ ∣ ∂f ∂y (t, y)
∣∣ ≤ L for all (t, y) ∈ R.
Liptischitz Condition and Well-Posednes
Definition. An IVP is said to be well-posed if a small perhubation in the data of the problem leads to only a small change in the solution.
Since numerical computation may very well introduce some perhubations to the problem, it is important that the problem that is to be solved is well-posed. Fortunately, the Lipschitz condition is a sufficient condition for the IVP problem to be well- posed.
Theorem (Well-Posedness of the IVP problem). If f (t, y) Satisfies the Lipschitz Condition, then the IVP is well-posed.
One of the simplest methods for solving the IVP is the classical Euler method.
The method is derived from the Taylor Series expansion of the function y(t). The function y(t) has the following Taylor series expansion of order n at t = ti+1:
y(ti+1) = y(ti) + (ti+1 − ti)y′(ti) + (ti+1 2!^ −^ ti)
2 y′′(ti) + · · · + (ti+1 n^ −!^ ti)
n y(n)(ti) + (ti+1 − ti) (n + 1)!
n+ yn+1(ξi), where ξi is in (ti, ti+1).
Substitute h = ti+1 − ti. Then Taylor Series Expansion of y(t) of order n at t = ti+1: y(ti+1) = y(ti) + hy′(ti) + h 2 2! y
′′(ti) + · · · + hn n! y
(n)(ti) + hn+ (n + 1)! y
(n+1)(ξi).
For n = 1, this formula reduces to
y(ti+1) = y(ti) + hy′(ti) + h
2 2 y
′′(ξ).
Algorithm: Euler’s Method for IVP Input: (i). The function f (t, y) (ii). The end points of the interval [a, b] : a and b (iii). The initial value: α = y(t 0 ) = y(a) Output: Approximations yi+1 of y(ti + 1), i = 0, 1 , · · · , N − 1. Step 1. Initialization: Set t 0 = a, y 0 = y(t 0 ) = y(a) = α. and N = b^ − h a. Step 2. For i = 0, 1 , · · · , N − 1 do Compute yi+1 = yi + hf (ti, yi) End
Example: y′^ = t^2 + 5, 0 ≤ t ≤ 1. y(0) = 0, h = 0.25. Find y 1 , y 2 , y 3 , and y 4 , approximations of y(0.25), y(0.50, y(0.75), and y(1), respectively. The points of subdivisions are: t 0 = 0, t 1 = 0. 25 , t 2 = 0. 50 , t 3 = 0.75 and t 4 = 1.
i = 0 : t 1 = t 0 + h = 0. 25 y 1 = y 0 + hf (t 0 , y 0 ) = 0 + .25(5) = 1. 25 (y(0.25) = 1.2552). i = 1 : t 2 = t 1 + h = 0. 50 y 2 = y 1 + hf (t 1 , y 1 ) = 1.25 + 0.25(t^21 + 5) = 1.25 + 0.25((0.25)^2 + 5) = 2. 5156 (y(0.5) = 2.5417). i = 2 : t 3 = t 2 + h = 0. 75 y 3 = y 2 + hf (t 2 , y 2 ) = 2.5156 + .25((.5)^2 + 5) = 3. 8281 (y(0.75) = 3.8906). etc.
Example: y′^ = t^2 + 5, 0 ≤ t ≤ 2 , y(0) = 0, h = 0. 5
So, the points of subdivisions are: t 0 = 0, t 1 = 0. 5 , t 2 = 1, t 3 = 1. 5 , t 4 = 2.
We compute y 1 , y 2 , y 3 , and y 4 , which are, respectively, approximations to y(0.5), y(1), y(1.5), and y(2).
after two terms. Thus, in obtaining Euler’s method, the first term neglected is h
2 2 y ′′(t).
So the local error in Euler’s method is: ELE = h
2 2 y
′′(ξi),
where ξi lies between ti and ti+1. In this case, we say that the local error is of order h^2 , written as O(h^2 ). Note that the local error ELE converges to zero as h → 0.
Global error is the difference between the true solution y(ti) and the approximate solution yi at t = ti. Thus, Global error = y(ti) − yi. Denote this by EGE.
The following theorem shows that the global error, EGE , is of order h.
Theorem: (Global Error Bound for the Euler Method) Let y(t) be the unique solution of the IVP: y′^ = f (t, y); y(a) = α. a ≤ t ≤ b, −∞ < y < ∞, Let L and M be two numbers such that ∣∣ ∣∣^ ∂f^ ∂y(t, y)
∣∣ ≤ L, and |y′′(t)| ≤ M in [a, b]. Then the global error EGE at t = ti satisfies |EGE | =
∣∣y(ti) − yi
∣∣ ≤ hM 2 L (e
L(ti−a) (^) − 1).
Thus, The global error bound for Euler’s method depends upon h, whereas the local error depends upon h^2.
Proof of the above theorem can be found in the book by G.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice Hall, Inc. (1971).
Remark. Since the exact solution y(t) of the IVP is not known, the above bound may not
be of practical importance as far as knowing how large the error can be a priori. However, from this error bound, we can say that the Euler method can be made to converge faster by decreasing the step-size. Furthermore, if the equalities, L and M of the above theorem can be found, then we can determine what step-size will be needed to achieve a certain accuracy, as the following example shows.
Example: dy dt = t
(^2) + y 2 2 , y(0) = 0 0 ≤ t ≤ 1 , − 1 ≤ y(t) ≤ 1. Determine how small the step-size should be so that the error does not exceed = 10−^4. Compute L: Since f (t, y) = t (^2) + y 2 2 ,^ we^ have ∂f ∂y =^ y Thus,
∣∣^ ∂f ∂y
∣∣ ≤ 1 for all y, giving L=. Find M:
Thus, y(ti+1) = y(ti)+hf (ti, y(ti))+ h
2 2 f^
′(ti, y(ti))+· · ·+ hn−^1 (n − 1)! f^
(n−1)(ti, yi)+ hn n! f^
(n−1)(ti, y(ti))+ hn+ (n + 1)! f^
(n−1)(ξi, y(ξi)
= y(ti) + h
f (ti, y(ti)) + h 2 f ′(ti, y(ti)) + · · · + h n− 1 n! f^ n− (^1) (ti, y(ti))
Tk(ti, y 1 ) = f (ti, yi) + h 2 f ′(ti, yi) + · · · + h k− 1 k! f^ (k−1)(ti, yi)
So, if we truncate the Taylor Series after (k + 1) terms and use the truncated series to obtain the approximating of y1+1 of y(ti+1), we have the following of k-th order Taylor’s algo- rithm for the IVP.
Taylor’s Algorithm of order k for IVP Input: (i) The function f (t, y) (ii) The end points: a and b (iii) The initial value: α = y(to) = y(a) (iv) The order of the algorithm: k (v) The step size: h Step 1 Initialization: t 0 = a, y 0 = α, N = b^ − h^ a Step 2. For i =... , N − 1 do 2.1 Compute Tk(ti, yi) = f (ti, yi) + h 2 f ′(ti, yi) +... + h
k− 1 k! f^
(k−1)(ti, yi) 2.2 Compute yi+1 = yi + hTk(ti, yi) End
Note: With k = 1, the above formula for yi+1, reduces to Euler’s method.
Example: y′^ = y − t^2 + 1, 0 ≤ t ≤ 2 , y(0) = 0. 5 , h = 0 · 2.
The points of division are: t 0 = 0, t 1 = 0.2, t 2 = 0.4, t 3 = 0.6, t 4 = 0.8, t 5 = 1, and so on.
Suppose that we want an expression of the approximation yi+1 in the form:
yi+1 = yi + α 1 k 1 + α 2 k 2 , (4.1) where k 1 = hf (ti, yi), (4.2)
and k 2 = hf (ti + αh, yi + βk 1 ). (4.3)
The constants α 1 and α 2 and α and β are to be chosen so that the formula is as accurate as the Taylor’s Series Method of as high as possible.
To develop the method we need an important result from Calculus: Taylor’s series for function to two variables.
Taylor’s Theorem for Function of Two Variables Let f (t, y) and its partial derivatives of orders up to (n + 1) are continuous in the domain D = {(t, y)|a ≤ t ≤ b, c ≤ y ≤ d}. Then f (t, y) = f (t 0 , y 0 ) +
(t − t 0 )∂f ∂t (t 0 , y 0 ) + (y − y 0 )∂f ∂y (t 0 , y 0 )
n!
∑^ n h=
(n i
(t − t 0 )h−i(y − y 0 )i^ ∂ nf ∂tn−^1 ∂yi^ (t^0 , y^0 )
Using the above theorem with n = 1, we have
f (ti + αh, yi + βk 1 ) = f (ti, yi) + αh ∂f ∂t (ti, yi) + βk 1 ∂f ∂y (ti, yi) (4.4)
From (4.4) and (4.3), we obtain
k 2 h =^ f^ (ti, yi) +^ αh
∂f ∂t (ti.yi) +^ βk^1
∂f ∂y (ti, yi).^ (4.5) Again, substituting the value of k 1 from (4.2) and k 2 from (4.3) in (4.1) we get (after some rearrangment):
yi+1 = yi + α 1 hf (ti, yi) + α 2 h
f (ti, yi) + αh ∂f ∂t (ti, yi) + βhf (ti, yi)∂f ∂y (ti, yi)
= yi + (α 1 + α 2 )hf (ti, yi) + α 2 h^2
α ∂f ∂t (ti, yi) + βf (ti, yi)∂f ∂y (ti, yi)
Also, note that y(ti+1) = y(ti) + hf (ti, yi) + h
2 2
(∂f ∂t (ti, yi) +^ f^ (ti, yi)
∂f ∂y (ti, yi)
So, neglecting the higher order terms, we can write
yi+1 = yi + hf (ti, yi) + h 2 2
(∂f ∂t (ti, yi) +^ f
∂f ∂y (ti, yi)
If we want (4.6) and (4.7) to agree for numerical approximations, then we must have
Since the number of unknowns here exceeds the number of equations, there are infinitely many possible solutions. The simplest solution is:
Algorithm: The Modified Euler Method Inputs: The given function: f (t, y) The end points of the interval: a and b The step-size: h The initial value y(t 0 ) = y(a) = α Outputs: Approximations yi+1 of y(ti+1) = y(t 0 + ih), i = 0, 1 , 2 , · · · , N − 1 Step 1 (Initialization) Set t 0 = a, y 0 = y(t 0 ) = y(a) = α N = b^ − h^ a Step 2 For i = 0, 1 , 2 , · · · , N − 1 do Compute k 1 = hf (ti, yi) Compute k 2 = hf (ti + h, yi + k 1 ) Compute yi+1 = yi +^12 (k 1 + k 2 ). End Example: y′^ = et, y(0) = 1, h = 0.5, 0 ≤ t ≤ 1 t 0 = 0, t 1 = 0.5, t 2 = 1 i = 0 : k 1 = hf (t 0 , y 0 ) = 0. 5 et^0 = 0. 5 k 2 = hf (t 0 + h, y 0 + k 1 ) = 0.5(et^0 +h) = 0. 5 e^0.^5 = 0. 8244 y 1 = y 0 + 12 (k 1 + k 2 ) = 1 + 0.5(0.5 + 0.8244) = 1. 6622 (y(0.5) = e^0.^5 = 1.6487) i = 1 : k 1 = hf (t 1 , y 1 ) = 0. 5 et^1 = 0. 5 e^0.^5 = 0. 8244 k 2 = hf (t 1 + h, y 1 + k 1 ) = 0. 5 e^0 .5+0.^5 = 0. 5 e = 1. 3591 y 2 = y 1 + 12 (k 1 + k 2 ) = 1.6622 + 12 (0.8244 + 1.3591) = 2. 7539 (y(1) = 2 .7183). Example: Given: y′^ = t + y, y(0) = 1, compute y 1 (approximation to y(0.01)) and y 2 (approximation to y(0.02) by using Modified Euler Method.
h = 0. 01 , y 0 = y(0) = 1. i = 0 : y 1 = y 0 +^12 (k 1 + k 2 ) k 1 = hf (t 0 , y 0 ) = 0.01(0 + 1) = 0. 01 k 2 = hf (t 0 + h, y 0 + k 1 ) = 0. 01 × f (0. 01 , 1 + 0.01) = 0. 01 × (0.01 + 1.01) = 0. 01 × 1 .02 = 0. 0102 Thus y 1 = 1 +^12 (0.01 + 0.0102) = 1. 0101 (Approximate value of y(0.01)
i = 1 : y 2 = y 1 +^12 (k 1 + k 2 ) k 1 = hf (t 1 , y 1 ) = 0. 01 × f (0. 01 , 1 .0101) = 0. 01 × (0.01 + 1.0101) = 0. 0102 k 2 = hf (t 1 + h, y 1 + k 1 ) = 0. 01 × f (0. 02 , 1 .0101 + 0.0102) = 0. 01 × (0.02 + 1.0203) = − 0. 0104
y 2 = 1.0101 +^12 (0.0102 + 0.0104) = 1.0204 (Approximate value of y(0.02)). Local Error in the Modified Euler Method
Since in deriving the modified Euler method, we neglected the terms involving h^3 and higher powers of h, the local error for this method is O(h^3 ). Thus with the modified Euler method, we will be able to use larger step-size h than the Euler method to obtain the same accuracy.
The Midpoint and Heun’s Methods
In deriving the modified Euler’s Method, we have considered only one set of possible values of α 1 , α 2 , α 1 and β. We will now consider two more sets of values.
