**Chapter 5**

**Initial-Value Problems for Ordinary** **Differential Equations**

**Hung-Yuan Fan (范洪源)**

**Department of Mathematics,**
**National Taiwan Normal University, Taiwan**

**Spring 2016**

**Section 5.1**

**The Elementary Theory of** **Initial-Value Problems**

**(初值問題的基本理論)**

**Objectives**

Develop numerical methods for approximating the solution to

**initial-value problem (IVP)**

{ _{dy}

*dt* *= f(t, y),* *a≤ t ≤ b,*

*y(a) = α,* (1)

**where y(t) is the unique solution to IVP (1) on [a, b].**

Error analysis for these numerical methods.

Note:

**1** **The first equation in (1) is an ordinary differential equation**

**(ODE; 常微分⽅程式).**

**2** **y(a) = α is called an initial condition (IC; 初值條件).**

**Def 5.1, p. 261**

**A function f(t, y) satisfies a Lipschitz condition in y on a set***D⊆ R*^{2} if**∃ a Lipschitz constant L > 0 s.t.**

*|f(t, y*1)*− f(t, y*2)| ≤ L|y1*− y*2*|,*
*whenever (t, y*_{1})*∈ D and (t, y*2)*∈ D.*

**Thm 5.4 (IVP 解的唯⼀性)**

*Suppose that f(t, y) is conti. on D ={(t, y) | a ≤ t ≤ b and y ∈ R}.*

If*f satisfies a Lipschitz condition in y on D, then the IVP (1) has a*

**unique solution y(t) for a**

**unique solution y(t) for a**

*≤ t ≤ b.*

**Corollary of Thm 5.4**

*Suppose that f(t, y) is conti. on D ={(t, y) | a ≤ t ≤ b and y ∈ R}.*

If*∃ a Lipschitz constant L > 0 with*
*∂f*

*∂y(t, y)** ≤ L* *∀ (t, y) ∈ D,*

then*f(t, y) satisfies a Lipschitz condition in y on D, and therefore,*
**the IVP (1) has a unique solution y(t) for a**≤ t ≤ b.

**Def 5.5, p. 263 (Well-Posedness of IVP)**

The IVP { _{dy}

*dt* *= f(t, y),* *a≤ t ≤ b,*
*y(a) = α*

**is said to be a well-posed problem if**

**A unique solution y(t) exists on [a, b], and**

*∃ ε*0 *and k > 0 s.t.for any 0 < ε < ε*_{0}*, whenever δ(t)∈ C[a, b]*

with *|δ(t)| < ε* *for t∈ [a, b] and* *|δ*0*| < ϵ***, the perturbed IVP**
{ _{dz}

*dt* *= f(t, z) +δ(t),* *a≤ t ≤ b,*
*z(a) = α +δ*0

**has a unique solution z(t) satisfying**

**Thm 5.6 (初值問題是 Well-Posed 的充分條件)**

*Suppose that D ={(t, y) | a ≤ t ≤ b and − ∞ < y < ∞} ⊆ R*^{2}.
If*f∈ C(D) satisfies a Lipschitz condition in y on D*, thenthe IVP
**(1) is well-posed.**

**Remarks**

Because anyround-off errorintroduced in the representation
perturbs the original IVP (1), numerical methods will always
**be connected with solving a perturbed IVP.**

**If the original IVP is well-posed, the numerical solution to a**
perturbed problem will accurately approximate the unique
solution to the original problem!

**Section 5.1 勾選習題**

**1, 3**

**Section 5.2**

**Euler’s Method**

**(尤拉法或歐拉法)**

### Derivation of Euler’s Method

Assume that the following IVP (1)
{ _{dy}

*dt* *= f(t, y),* *a≤ t ≤ b,*
*y(a) = α*

**is well-posed and y(t) is the unique sol. to IVP (1) on [a, b].**

*Choose the equally-distributed mesh points (網格點) on [a, b]*

*t*_{i}*= a + i· h, i = 0, 1, 2, . . . ,N,* (2)
where*N∈ N*and *h = (b− a)/N* **is the step size. (步⻑)**
**Note**

*The graph of the unique solution y(t) evaluated at each mesh*
point is shown below.

### Derivation of Euler’s Method (Conti’d)

If*y(t)∈ C*^{2}*[a, b], it follows from Taylor’s Thm that for each*
*i = 0, 1, . . . , N− 1, ∃ ξ**i**∈ (t**i**, t**i+1*) s.t.

*y(t*_{i+1}*) = y(t**i**) + (t**i+1**− t**i**)y*^{′}*(t**i*) + *(t**i+1**− t**i*)^{2}
2 *y*^{′′}*(ξ**i*)

*= y(t**i*) +*hy*^{′}*(t**i*) +*h*^{2}
2 *y*^{′′}*(ξ**i*)

*= y(t** _{i}*) +

*hf(t*

_{i}*, y(t*

*)) +*

_{i}*h*

^{2}2

*y*

^{′′}*(ξ*

*) where*

_{i}*h*

*is the step size and t*

*is chosen as in (2).*

_{i}Deleting the remainder term =**⇒ Euler’s method constructs***w*_{i}*≈ y(t**i**) for i = 0, 1, . . . , N*by

The fist step of Euler’s method is shown below.

*After N steps of Euler’s method defined as in (3), the differences*
between*y(t** _{i}*)

*and w*

_{i}*(i = 1, 2, . . . , N) are shown below.*

**Algorithm 5.1: Euler’s Method**

INPUT *endpoints a, b; positive integer N; initial condition α.*

OUTPUT approximation *w* *to y at the (N + 1) values oft.*

Step 1 *Set h = (b− a)/N; t = a; w = α;*

*OUTPUT(t, w).*

Step 2 **For i = 1, 2, . . . , N do Steps 3–4.**

Step 3 Set*w = w + h**· f(t, w)**; (Compute w**i*.)
*t = a + i**· h**. (Compute t**i*.)
Step 4 *OUTPUT(t, w).*

Step 5

**STOP.**

**Example 1, p. 268**
Consider the following IVP

*y*^{′}*= y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

(a) Show that *y(t) = (t + 1)*^{2}*− 0.5e** ^{t}* is the unique solution to
above IVP.

(b) Apply Euler’s method (Alg. 5.1) with*h = 0.2* and*N = 10* to
obtain approximations *w** _{i}*, and compare these with the actual

*values of y(t*

_{i}*) for i = 1, 2, . . . , N.*

**Solution**

(a) *Note that y(0) = (1 + 0)*^{2}*− 0.5e*^{0} *= 0.5 and it is easily seen*
*that y(t) satisfies the given ODE by direct computation!*

(b) From the Euler’s method in (3), we have
*w*_{0} *= 0.5,* *t*_{0}*= 0,*

*w*_{i+1}*= w**i**+ 0.2(w**i**+ t*^{2}_{i}*− 1) = w**i**+ 0.2[w**i**− (0.2i)*^{2}+ 1]

=**1.2w**_{i}**− 0.008i**^{2}**+ 0.2,***i = 0, 1, . . . , 9,*
where*t*_{i}*= 0 + 0.2i = 0.2ifor i = 0, 1, 2, . . . , 10 = N.*

The numerical results of Part (b) are shown in the following table.

### Error Bounds for Euler’s Method

**Thm 5.9 (Theoretical Error Bound)**

Suppose that*f is conti.* and satisfiesa Lipschitz condition with
*constant L > 0on D ={(t, y) ∈ R*^{2}*| a ≤ t ≤ b, −∞ < y < ∞},*
and that*∃ M > 0* with

*|y*^{′′}*(t)| ≤ M ∀ t ∈ [a, b],*

*where y(t) is the unique sol. to the IVP (1). If w*_{0}*, w*_{1}*,· · · , w**N* are
*the approximations generated by Euler’s method for some N∈ N,*
*then for each i = 0, 1, . . . , N, we have*

*|y(t**i*)*− w**i**| ≤* *hM*
*2L*
[

*e*^{L(t}^{i}^{−a)}*− 1*]
*with t*_{i}*being the grid points and h being the step size.*

**Comments on Thm 5.9**

*In practice, it is difficult to verify the boundedness of y*^{′′}*(t)!*

Instead, we may check for the boundednes of
*y*^{′′}*(t) =* *dy*^{′}*(t)*

*dt* = *d*

*dtf(t, y(t))*

= *∂f*

*∂t(t, y(t)) +* *∂f*

*∂y(t, y(t))· f(t, y(t))*
**without explicitly knowing the unique solution y(t).**

### A Test for Theoretical Error Bound

**Example 2, p. 272 (驗證 Thm 5.9 的誤差上界)**

**As in Example 1, Euler’s method with***h = 0.2*is applied for
computing the approximations*w*_{i}*(i = 0, 1, . . . , N)* of the unique
solution*y(t) = (t + 1)*^{2}*− 0.5e** ^{t}* to the IVP

*y*^{′}*= y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

Compare the error bounds given in Thm 5.9 to the actual errors

*|y(t**i*)*− w**i**| for i = 0, 1, . . . , N.*

**Solution (1/2)**

*Let f(t, y) = y− t*^{2}+ 1 be a real-valued function defined on the set
*D ={(t, y) ∈ R*^{2}*| 0 ≤ t ≤ 2, −∞ < y < ∞}.*

Then*f∈ C(D) satisfies a Lipschitz condition in y on D with L = 1*,
since

*∂f*

*∂y(t, y) = 1* or *∂f*

*∂y(t, y)** ≤ 1 ∀ (t, y) ∈ D.*

*Moreover, since the unique sol. is y(t) = (t + 1)*^{2}*− 0.5e** ^{t}*, we have

*y*

^{′′}*(t) = 2− 0.5e*

*and hence*

^{t}*|y*^{′′}*(t)| ≤0.5e*^{2}*− 2 ≡ M* *∀ t ∈ [0, 2].*

**Solution (2/2)**

So, it follows from Thm 5.9 that the error bounds for Euler’s method are given by

*|y(t**i*)*− w**i**| ≤* *hM*
*2L*

[

*e*^{L(t}^{i}^{−a)}*− 1*]

*= (0.1)· (0.5e*^{2}*− 2) · (e*^{t}^{i}*− 1),*
*where the approx. w** _{i}* computed by Euler’s method are

*w*_{0} *= 0.5,* *w*_{i+1}*= 1.2w**i**− 0.008i*^{2}*+ 0.2 for i = 0, 1, . . . , 9,*
*and the mesh points are t*_{i}*= 0.2i for i = 0, 1, 2, . . . , 10 = N.*

The numerical comparison between actual errors and error bounds is shown in the following table.

**Rate of Convergence for Euler’s Method**

*|y(t**i*)*− w**i**| =O(h)* *for each i = 0, 1, . . . , N.*

**Finite-Digit Approximations to y(t*** _{i}*)

*If h = b − aN* *and t*_{i}*= a + ih for i = 0, 1, . . . , N, then note that*
Euler’s method is performed

**in the exact arithmetic:**

*w*_{0}*= α,*

*w*_{i+1}*= w**i**+ hf(t**i**, w**i*) *for i = 0, 1, . . . , N− 1.*

**in the finite-digit arithmetic:**

*u*_{0}*= α +δ*_{0}*,*

*u*_{i+1}*= u**i**+ hf(t**i**, u**i*) +*δ**i+1* *for i = 0, 1, . . . , N− 1,* (4)
where*δ**i* *denotes the round-off error associated with u** _{i}* for

*each i = 0, 1, . . . , N.*

### Practical Error Bounds for Euler’s Method

**Thm 5.10 (Error Bound in Finite-Digit Arithmetic)**

*Let y(t) be the unique solution to the IVP (1) and u*_{0}*, u*_{1}*,· · · , u**N*

be finite-digit approximations obtained using (4). If*|δ*_{i}*| < δ for*
*each i = 0, 1, . . . , N and the sufficient conditions of Thm 5.9 hold,*
then

*|y(t**i*)*− u**i**| ≤* 1
*L*

(*hM*
2 +*δ*

*h*
)[

*e*^{L(t}^{i}^{−a)}*− 1*]

+*|δ*0*|e*^{L(t}^{i}* ^{−a)}* (5)

*for each i = 0, 1, . . . , N.*

**Comments on Thm 5.10**
Since it is easily seen that

*h*lim*→0*^{+}

(*hM*
2 +*δ*

*h*
)

=*∞,*

**deceasing**

*h*

**tends to increase the total error in the** **approximation.**

If we let

*E(h) =* *hM*
2 + *δ*

*h* for *h> 0,*

*then E** ^{′}*(h) = M

*2 −*

^{h}

^{δ}^{2}, and therefore, it follows from Calculus

*that E(h) is minimized at*

*h*

*=*

^{∗}√*2δ*

*M*. In fact, we know that
*E*^{′}*(h) < 0* or *E(h) is decreasing for 0 < h < h*^{∗}*,*
*E*^{′}*(h) > 0* or *E(h) is increasing for h > h*^{∗}*.*

**Section 5.2 勾選習題**

**1, 3**

**Section 5.3**

**High er-Order Taylor Methods**

**(⾼階泰勒法)**

*Let y(t) be the unique solution to the IVP*

*y*^{′}*= f(t, y),* *a≤ t ≤ b, y(a) = α.* (6)

**Def 5.11 (局部截斷誤差)**

The difference method (差分⽅法) for solving the IVP (6)
*w*_{0}*= α,* *w*_{i+1}*= w**i**+ hϕ(t**i**, w**i*) *for i = 0, 1, . . . , N− 1,*
**has local truncation error (簡稱 LTE)**

*τ**i+1**(h) =* *y*_{i+1}*− (y**i**+ hϕ(t**i**, y**i*))

*h* = *y*_{i+1}*− y**i*

*h* *− ϕ(t**i**, y**i*),

**The LET of Euler’s Method**

If*y(t)∈ C*^{2}*[a, b], it follows from Taylor’s Thm that for each*
*i = 0, 1, . . . , N− 1, ∃ ξ**i**∈ (t**i**, t**i+1*) s.t.

*y*_{i+1}*= y**i**+ hf(t**i**, y**i*) +*h*^{2}

2 *y*^{′′}*(ξ**i**),* (7)
*where h is the step size and t** _{i}* is chosen as in (2).

*From (7), LTE of Euler’s method at the ith step is*
*τ**i+1**(h) =* *y*_{i+1}*− y**i*

*h* *−ϕ(t**i**, y**i*) = *y*_{i+1}*− y**i*

*h* *−f(t**i**, y**i*) = *h*
2*y*^{′′}*(ξ**i**).*

So, we see that

*τ**i**(h) = O(h)* *for each i = 1, 2, . . . , N,*
since*y*^{′′}*is bounded on [a, b].*

*Taylor Method of Order n* *∈ N*

*If the sol. y(t) is smooth enough, say,* *y∈ C*^{n+1}*[a, b], then*

*∃ ξ**i**∈ (t**i**, t**i+1*) s.t.

*y(t*_{i+1}*) = y(t**i*) +

∑*n*
*k=1*

*h*^{k}

*k!y*^{(k)}*(t**i*) + *h*^{n+1}

*(n + 1)!y*^{(n+1)}*(ξ**i*) (8)
*for each i = 0, 1, . . . , N− 1.*

Successive differentiation gives that

*y*^{′}*(t) = f(t, y(t)), y*^{′′}*(t) = f*^{′}*(t, y(t)),· · · , y*^{(n+1)}*(t) = f*^{(n)}*(t, y(t)).*

Then Eq. (8) can be rewritten as

∑*n* *h*^{k}^{−1}_{−1)}*h*^{n+1}

**Taylor’s Method of Order n**∈ N

*The approximations w*_{i}*to y(t*_{i}*) (i = 0, 1, . . . , N) are computed by*
*w*_{0}*= α,*

*w*_{i+1}*= w**i**+ hT*^{(n)}*(t**i**, w**i*) *for each i = 0, 1, . . . , N− 1,* (10)
*where h is the step size and*

*T*^{(n)}*(t*_{i}*, w*_{i}*) = f(t*_{i}*, w** _{i}*) +

*h*

2*f*^{′}*(t*_{i}*, w** _{i}*) +

*· · · +*

*h*

^{n}

^{−1}*n!* *f*^{(n}^{−1)}*(t*_{i}*, w** _{i}*)

*for each i = 0, 1, . . . , N− 1.*

**Note:**

Euler’s method is just the Taylor’s method of order one!
**Thm 5.12 (⾼階泰勒法的局部截斷誤差)**

*Let y(t) be the unique solution to the IVP (6) on [a, b]. If*
*y∈ C*^{n+1}*[a, b], then the LTEs of Taylor’s method of order n*
defined in (10) satisfy

*τ**i**(h) = O(h** ^{n}*)

*for each i = 1, 2, . . . , N,*

*where n and N are some positive integers.*

**Recall from Eq. (9)**

*For each i = 0, 1, . . . , N− 1, ∃ ξ**i**∈ (t**i**, t**i+1*) s.t.

*y* *·*

∑*n* *h*^{k}^{−1}

*f*^{(k}^{−1)}*h*^{n+1}*f*^{(n)}

**Proof**

From Taylor’s Thm and (9), we obtain
*y*_{i+1}*− y**i*

*h* *−T*^{(n)}*(t**i**, y**i*)= *h*^{n}

*(n + 1)!f*^{(n)}*(ξ**i**, y(ξ**i*))
*for some ξ*_{i}*∈ (t**i**, t**i+1*).

Sine*y*^{(n+1)}*= f*^{(n)}*∈ C[a, b]*,*|f*^{(n)}*(t, y(t))| < M ∀ t ∈ [a, b]*

*and hence the LTE at the ith step satisfies*

*|τ**i+1**(h)| ≤* *M*

*(n + 1)!h*^{n}*for h > 0,*
i.e., *τ**i+1**(h) = O(h*^{n}*) for each i = 0, 1, . . . , N− 1.*

**Example 1, p. 278**

**Apply Taylor’s method of orders (a) two and (b) four with**
*N = 10*to compute the approximations*w*_{i}*(i = 0, 1, . . . , N)*of the
unique solution

*y(t) = (t + 1)*^{2}*− 0.5e** ^{t}*
to the IVP

*y*^{′}*= y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

*Let f be a real-valued function defined by*

*f(t, yt)) = y(t)− t*^{2}+ 1 *∀ t ∈ [0, 2].*

**Solution of Part (a): Taylor’s Method of Order 2**
*Since f*^{′}*(t, y) = y*^{′}*− 2t = y − t*^{2}+ 1*− 2t, we have*

*T*^{(2)}*(t**i**, w**i**) = f(t**i**, w**i*) +*h*

2*f*^{′}*(t**i**, w**i*)

= (1 +*h*

2*)(w**i**− t*^{2}*i* + 1)*− ht**i*

*for each i = 0, 1, . . . , N− 1.*

Taylor’s method of order 2 is defined by
*w*_{0}*= 0.5,*

*w*_{i+1}*= w**i**+ hT*^{(2)}*(t**i**, w**i*)

*= w*_{i}*+ (0.2)*
[

(1 +*0.2*

2 *)(w*_{i}*− (0.2i)*^{2}+ 1)*− (0.2)(0.2i)*]

**= 1.22w***i***− 0.0088i**^{2}**− 0.008i + 0.22***for each i = 0, 1, . . . , 9.*

The numerical results of Part (a) are shown in the following table.

**Solution of Part (b): Taylor’s Method of Order 4 (1/2)**
*Similarly, successive differentiation w.r.t. t gives*

*f*^{′′}*(t, y) = y− t*^{2}*− 2t − 1 = f*^{(3)}*(t, y).*

*The function T*^{(4)}*(t**i**, w**i*) is defined by
*T*^{(4)}*(t**i**, w**i**) = f(t**i**, w**i*) +*h*

2*f*^{′}*(t**i**, w**i*) +*h*^{2}

6*f*^{′′}*(t**i**, w**i*)
+*h*^{3}

24*f*^{(3)}*(t**i**, w**i*)

= (1 +*h*
2 +*h*^{2}

6 + *h*^{3}

24*)(w**i**− t*^{2}*i*)*− (1 +* *h*
3 +*h*^{2}

12*)(ht**i*)
+ 1 +*h*

2 *−h*^{2}
6 *−* *h*^{3}

24*.*

**Solution of Part (b): Taylor’s Method of Order 4 (2/2)**
*Substituting h = 0.2 and t*_{i}*= 0.2i into T*^{(4)}, we thus obtain
Taylor’s method of order 4 as

*w*_{0}*= 0.5,*

*w*_{i+1}*= w**i**+ hT*^{(4)}*(t**i**, w**i*)

**= 1.2214w**_{i}**− 0.008856i**^{2}**− 0.00856i + 0.2186***for each i = 0, 1, . . . , 9.*

The numerical results of Part (b) are shown in the following table.

**Section 5.3 勾選習題**

**5, 7**

**Section 5.4**

**Runge-Kutta Methods**

**(R-K 法)**

**Taylor’s Methods v.s. Runge-Kutta Methods**

**1** Taylor’s Methods

**Advantage: They have hogh-order LTE at each step.**

**Disadvantage: They require computation and evaluation of**
*the high-order derivatives of f(t, y(t)) w.r.t. the variable t.*

**2** Runge-Kutta Methods

They also have thehigh-order LTE at each step.

But eliminate the need to compute and evaluate the
*derivatives of f(t, y).*

**Thm 5.13 (Taylor’s Thm for Functions of Two Variables)**
*Suppose f(t, y) and all its partial derivatives of order≤ n + 1 are*
*conti. on D ={(t, y) | a ≤ t ≤ b, c ≤ y ≤ d}, and let*

*Q*_{0}*(t*_{0}*, y*_{0})*∈ D. Then for every (t, y) ∈ D, ∃ ξ between t and t*0 and

*∃ µ between y and y*0 s.t.

*f(t, y) = P*_{n}*(t, y) + R**n**(t, y),*
where

*P*_{n}*(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)
]

+

[*(t− t*0)^{2}

2 *f*_{tt}*(Q*_{0}*) + (t− t*0*)(y− y*0*)f**yt**(Q*_{0}) +*(y− y*0)^{2}

2 *f*_{yy}*(Q*_{0})
]

+*· · · +*[1
*n!*

∑*n*
*j=0*

( *n*
*j*

)

*(t− t*0)^{n}^{−j}*(y− y*0)^{j}*∂*^{n}*f*

*∂t*^{n}^{−j}*∂y*^{j}*(t*0*, y*0)
]

(11)

**Thm 5.13–Conti’d**
and

*R**n**(t, y) =* 1
*(n + 1)!*

∑*n+1*
*j=0*

( *n + 1*
*j*

)

*(t**−t*0)^{n+1}^{−j}*(y**−y*0)^{j}*∂*^{n+1}*f*

*∂t*^{n+1}^{−j}*∂y*^{j}*(ξ, µ).*

(12)
*In this case, P*_{n}**(t, y) is called the nth Taylor polynomial in two**

**variables for the function f about Q**

**variables for the function f about Q**

_{0}

*(t*0

*, y*0

*), and R*

*n*

*(t, y) is the*

**remainder term associated with P**

**remainder term associated with P**

_{n}*(t, y).*

**A Class of R-K Methods**

**1** Runge-Kutta Methods of Order 2 (⼆階 R-K 法)
Midpoint method (中點法)

Modified Euler method (修正歐拉法)
LTE at each step =*O(h*^{2})

**2** Runge-Kutta Methods of Order 3 (三階 R-K 法)
Heun’s method

LTE at each step =*O(h*^{3})

**3** Runge-Kutta Methods of Order 4 with LTE*O(h*^{4})

### The Midpoint Method (1/2)

**For the Midpoint Method, we try to determine the values**
*for a*_{1}*, α*_{1} *and β*_{1} s.t. *a*_{1}*f(t + α*_{1}*, y + β*_{1}) approximates

*T*^{(2)}*(t, y) = f(t, y) +h*
2*f*^{′}*(t, y)*

*= f(t, y) +h*
2

[*f*_{t}*(t, y) + f**y**(t, y)· y*^{′}*(t)*
]

*= f(t, y) +h*

2*f*_{t}*(t, y) +h*

2*f*_{y}*(t, y)· f(t, y)* (13)
From Thm 5.13 with (11)–(12) *⇒ ∃ ξ between t + α*1 *and t*
and*∃ µ between y + β*1 *and y s.t.*

*a*_{1}*f(t + α*_{1}*, y + β*1)*= a*1*f(t, y)+(a*_{1}*α*1*)f**t**(t, y)+(a*1*β*1*)f**y**(t, y)+R*1*,*
(14)

### The Midpoint Method (2/2)

Comparing the coefficients in (13) and (14) *⇒*
*a*_{1} *= 1,* *α*1 = *h*

2*,* *β*1 = *h*
2*f(t, y).*

So, we further obtain

*T*^{(2)}*(t, y) = f(t +h*
2*, y +* *h*

2*f(t, y))− R*1

*with R*_{1}= ^{h}_{8}^{2}*f*_{tt}*(ξ, µ) +* ^{h}_{4}^{2}*f(t, y)f*_{yt}*(ξ, µ) +*^{h}_{8}^{2}*[f(t, y)]*^{2}*f*_{yy}*(ξ, µ).*

**This means that the ith LTE of Midpoint Method satisfies***τ**i+1**(h) =* *y*_{i+1}*− y**i*

*h* *− f(t**i*+*h*

2*, y**i*+ *h*

2*f(t*_{i}*, y**i*))

= *y*_{i+1}*− y**i*

*h* *− T*^{(2)}*(t**i**, y**i**) + R*1=*O(h*^{2}),
since *R*_{1} *= O(h*^{2}*) if all second-order partial derivatives of f are*
*bounded on [a, b].*

**R-K Methods of Order Two**

**1**

**Midpoint Method:**

*w*_{0} *= α,*
*w*_{i+1}*= w**i**+ hf*

(*t** _{i}*+

*h*

2*, w**i*+*h*

2*f(t*_{i}*, w**i*)
)

*,*
*for each i = 0, 1, . . . , N− 1.*

**2**

**Modified Euler Method:**

*w*_{0} *= α,*
*w*_{i+1}*= w**i*+*h*

2 [

*f(t*_{i}*, w**i**) + f(t*_{i+1}*, w**i*+*hf(t*_{i}*, w**i*))
]

*,*

**Comments on R-K Methods of Order 2**

**Two function evaluations of f are required per step.**

**Two function evaluations of f are required per step.**

**For the Modified Euler Method, we use the function**
*a*_{1}*f(t, y) +a*_{2}*f(t +α*2*, y +δ*2*f(t, y))*

to approximate *T*^{(3)}*(t, y). The desired parameters are given by*
*a*_{1}*= a*_{2} = 1

2*,* *α*_{2} *= δ*_{2} *= h.*

The LTE at each step is *O(h*^{2}).

**Example 2, p. 286**

**Apply (a) Midpoint Method and (b) Modified Euler Method**
with*N = 10* to compute the approximations*w*_{i}*(i = 0, 1, . . . , N)*of
the unique solution

*y(t) = (t + 1)*^{2}*− 0.5e** ^{t}*
to the IVP

*y*^{′}*= f(t, y) = y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

The numerical results of Parts (a) and (b) are shown in the following table.

### R-K Methods of Order Three

The function*T*^{(3)}*(t, y)* can be approximated by
*f*(

*t +α*1*, y +δ*1*f*(

*t +α*2*, y +δ*2*f(t, y)*))
*,*
involving 4 parameters to be determined.

**Heun’s Method (R-K Method of Order 3)**

*w*_{0} *=α,*
*w*_{i+1}*=w**i*+*h*

4

[*f(t*_{i}*, w**i*)+

3f(
*t** _{i}*+

*2h*

*, w**i*+*2h*
*f*(

*t** _{i}*+

*h*

*, w**i*+*h*

*f(t*_{i}*, w**i*)
))]

*,*

### R-K Methods of Order Three

**Algorithm: Heun’s Method (Order Three)**

*w*_{0}*= α,*

*k*_{1}*= h· f(t**i**, w**i**),*
*k*_{2}*= h· f*(

*t** _{i}*+

*h*

3*, w**i*+1
3*k*_{1})

*,*
*k*_{3}*= h· f*(

*t** _{i}*+

*2h*

3 *, w**i*+2
3*k*_{2})

*,*
*w*_{i+1}*= w**i*+1

4*(k*_{1}*+ 3k*_{3}),
*for each i = 0, 1, . . . , N− 1.*

**Comments on R-K Methods of Order 3**

**Three function evaluations of f are required per step.**

**Three function evaluations of f are required per step.**

The LTE at each step is *O(h*^{3}).

R-K methods of order 3 are NOT generally used in practice!

**Example for Heun’s Method**

**Apply Heun’s Method with***N = 10*to compute the
approximations*w*_{i}*(i = 0, 1, . . . , N)* of the unique solution

*y(t) = (t + 1)*^{2}*− 0.5e** ^{t}*
to the IVP

*y*^{′}*= f(t, y) = y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

*The numerical results of Heun’s method with N = 10 and h = 0.2*
are shown in the following table.

### R-K Methods of Order Four

**Algorithm 5.2: Runge-Kutta (Order Four)**

*w*_{0} *= α,*

*k*_{1} *= h· f(t**i**, w**i**),*
*k*_{2} *= h· f*(

*t** _{i}*+

*h*

2*, w**i*+1
2*k*_{1})

*,*
*k*_{3} *= h· f*(

*t** _{i}*+

*h*

2*, w**i*+1
2*k*_{2})

*,*
*k*_{4} *= h· f(t*_{i+1}*, w**i*+*k*_{3}*),*
*w*_{i+1}*= w**i*+1

6*(k*1*+ 2k*2*+ 2k*3*+ k*4),
*for each i = 0, 1, . . . , N− 1.*

**Comments on R-K Methods of Order 4**

**Four function evaluations of f are required per step.**

**Four function evaluations of f are required per step.**

The LTE at each step is *O(h*^{4}).

These methods are the most commonly used for solving the IVPs in practice!

**Example 3, p. 289**

**Apply Runge-Kutta Method of Order 4 with***N = 10*to compute
the approximations*w*_{i}*(i = 0, 1, . . . , N)* of the unique solution

*y(t) = (t + 1)*^{2}*− 0.5e** ^{t}*
to the IVP

*y*^{′}*= f(t, y) = y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5.*

Numerical results for the R-K method of order 4 are shown in the following table.

**For R-K methods, the relationship between the number of**

**(function) evaluations per step and the order of LTE is shown**

in the following table.
**Reference**

*J. C. Butcher, The non-existence of ten-stage eight-order explicit*
**Runge-Kutta methods, BIT, Vol. 25, pp. 521–542, 1985.**

**Comparisons Between R-K Methods**

If the R-K method of order 4 is to be superior to Euler’s
method, it should give more accuracy*with step size h* than
Euler’s method *with step size h/4.*

If the R-K method of order 4 is to be superior to the R-K
method of order 2, it should give more accuracy with step size
*h* than second-order method *with step size h/2.*

### Illustration

For the same IVP as before

*y*^{′}*= f(t, y) = y− t*^{2}*+ 1,* 0*≤ t ≤ 2, y(0) = 0.5,*
**we apply Euler’s method with***h = 0.1/4 = 0.025, Modified*

**Euler method with**

*h = 0.1/2 = 0.05, and R-K method of order*

**4 with**

*h = 0.1. The numerical results are given as follows.*

**Section 5.4 勾選習題**