Numerical Solution of IVP for ODE

Introduction to the Numerical Solution of IVP for ODE

Consider the IVP: DE x

′

= f (t, x), IC x(a) = x

. For simplicity, we will assume here that

x(t) ∈ R

(so F = R), and that f (t, x) is continuous in t, x and uniformly Lipschitz in x

(with Lipschitz constant L) on [a, b] × R

. So we have global existence and uniqueness for

the IVP on [a, b].

Moreover, the solution of the IVP x

′

= f (t, x), x(a) = x

depends continuously on the

initial values x

∈ R

. This IVP is an example of a well-posed problem: for each choice of

the “data” (here, the initial values x

), we have:

(1) Existence. There exists a solution of the IVP on [a, b].

(2) Uniqueness. The solution, for each given x

, is unique.

(3) Continuous Dependence. The solution depends continuously on the data.

Here, e.g., the map x

7→ x(t, x

) is continuous from R

into (C ([a, b]) , k · k

∞

). A well-posed

problem is a reasonable problem to approximate numerically.

Grid Functions

· · · t

z}|{z}|{

Choose a mesh width h (with 0 < h ≤ b − a), and let
N =

b−a

(greatest integer ≤ (b−a)/h). Let t

= a + ih

(i = 0, 1, . . . , N) be the grid points in t (note: t

= a),

and let x

denote the approximation to x(t

). Note that

and x

depend on h, but we will usually suppress this

dependence in our notation.

Explicit One-Step Methods

Form of method: start with x

(presumably x

≈ x

). Recursively compute x

, . . . , x

i+1

= x

+ hψ(h, t

, x

i = 0, . . . , N − 1.

Here, ψ(h, t, x) is a function defined for 0 ≤ h ≤ b−a, a ≤ t ≤ b, x ∈ R

, and ψ is associated

with the given function f (t, x).

Examples.

i+1

solution through

(ti,x9)

tangent line

i+1

Euler’s Method.

i+1

= x

+ hf (t

, x

)

Here, ψ(h, t, x) = f (t, x).

Ordinary Differential Equations

Taylor Methods.

Let p be an integer ≥ 1. To see how the Taylor Method of order p is constructed, consider
the Taylor expansion of a C

p+1

solution x(t) of x

′

= f (t, x):

x(t + h) = x(t) + hx

′

(t) + · · · +

(p)

(t) +

O(h

p+1

)

| {z }

remainder term

where the remainder term is O(h

p+1

) by Taylor’s Theorem with remainder. In the ap-

proximation, we will neglect the remainder term, and use the DE x

′

= f (t, x) to replace

′

(t), x

′′

(t), . . . by expressions involving f and its derivatives:

′

(t) = f (t, x(t))

′′

(t) =

(f (t, x(t))) =

(n×1)

(t,x(t))

(n×n)

(t,x(t))

(n×1)

= (D

f + (D

f )f )

(t,x(t))

(for n = 1, this is f

+ f

f ).

For higher derivatives, inductively differentiate the expression for the previous derivative,
and replace any occurrence of

by f (t, x(t)). These expansions lead us to define the Taylor

methods of order p:

p = 1 :

i+1

= x

+ hf (t

, x

)

(Euler’s method, ψ(h, t, x) = f (t, x))

p = 2 :

i+1

= x

+ hf (t

, x

) +

f + (D

f )f )

)

For the case p = 2, we have

ψ(h, t, x) = T

(h, t, x) ≡

f +

f + (D

f )f )

(t,x)

We will use the notation T

(h, t, x) to denote the ψ(h, t, x) function for the Taylor method

of order p.

Remark.

Taylor methods of order ≥ 2 are rarely used computationally. They require deriva-

tives of f to be programmed and evaluated. They are, however, of theoretical interest in
determining the order of a method.

Remark.

A “one-step method” is actually an association of a function ψ(h, t, x) (defined for

0 ≤ h ≤ b − a, a ≤ t ≤ b, x ∈ R

) to each function f (t, x) (which is continuous in t, x and

Lipschitz in x on [a, b]×R

). We study “methods” looking at one function f at a time. Many

methods (e.g., Taylor methods of order p ≥ 2) require more smoothness of f, either for their

Introduction to the Numerical Solution of IVP for ODE

definition, or to guarantee that the solution x(t) is sufficiently smooth. Recall that if f ∈ C

(in t and x), then the solution x(t) of the IVP x

′

= f (t, x), x(a) = x

is in C

p+1

([a, b]).

For “higher-order” methods, this smoothness is essential in getting the error to be higher
order in h. We will assume from here on (usually tacitly) that f is sufficiently smooth when
needed.

Examples.

Modified Euler’s Method

i+1

= x

+ hf t

, x

f (t

, x

)

(so ψ(h, t, x) = f t +

, x +

f (t, x)

Here ψ(h, t, x) tries to approximate

′

t +

= f t +

, x t +

using the Euler approximation to x t +

≈ x(t) +

f (t, x(t))

Improved Euler’s Method

(or Heun’s Method)

i+1

= x

(f (t

, x

) + f (t

i+1

, x

+ hf (t

, x

)))

(so ψ(h, t, x) =

1
2

(f (t, x) + f (t + h, x + hf (t, x)))).

Here again ψ(h, t, x) tries to approximate

′

t +

≈

1
2

′

(t) + x

′

(t + h)).

Or ψ(h, t, x) can be viewed as an approximation to the trapezoid rule applied to

(x(t + h) − x(t)) =

t+h

′

≈

1
2

′

(t) +

1
2

′

(t + h).

Modified Euler and Improved Euler are examples of 2

order two-stage Runge-Kutta

methods. Notice that no derivatives of f need be evaluated, but f needs to be evaluated
twice

in each step (from x

to x

i+1

Before stating the convergence theorem, we introduce the concept of accuracy.

Local Truncation Error

Let x

i+1

= x

+ hψ(h, t

, x

) be a one-step method, and let x(t) be a solution of the DE

′

= f (t, x). The local truncation error (LTE) for x(t) is defined to be

l(h, t) ≡ x(t + h) − (x(t) + hψ(h, t, x(t))),

Ordinary Differential Equations

that is, the local truncation error is the amount by which the true solution of the DE fails to
satisfy the numerical scheme.

l(h, t) is defined for those (h, t) for which 0 < h ≤ b − a and

a ≤ t ≤ b − h.

Define

τ (h, t) =

l(h, t)

and set τ

(h) = τ (h, t

). Also, set

τ (h) = max

a≤t≤b−h

|τ(h, t)| for 0 < h ≤ b − a.

Note that

l(h, t

) = x(t

i+1

) − (x(t

) + hψ(h, t

, x(t

))),

explicitly showing the dependence of l on h, t

, and x(t).

Definition. A one-step method is called [formally] accurate of order p (for a positive integer
p) if for any solution x(t) of the DE x

′

= f (t, x) which is C

p+1

, we have l(h, t) = O(h

p+1

Definition. A one-step method is called consistent if ψ(0, t, x) = f (t, x). Consistency is
essentially minimal accuracy:

Proposition. A one-step method

i+1

= x

+ hψ(h, t

, x

where ψ(h, t, x) is continuous for 0 ≤ h ≤ h

, a ≤ t ≤ b, x ∈ R

for some h

∈ (0, b − a], is

consistent with the DE x

′

= f (t, x) if and only if τ (h) → 0 as h → 0

Proof. Suppose the method is consistent. Fix a solution x(t). For 0 < h ≤ h

, let

Z(h) =

max

a≤s,t≤b, |s−t|≤h

|ψ(0, s, x(s)) − ψ(h, t, x(t))|.

By uniform continuity, Z(h) → 0 as h → 0

. Now

l(h, t) = x(t + h) − x(t) − hψ(h, t, x(t))

t+h

′

(s) − ψ(h, t, x(t))] ds

t+h

[f (s, x(s)) − ψ(h, t, x(t))] ds

t+h

[ψ(0, s, x(s)) − ψ(h, t, x(t))] ds,

so |l(h, t)| ≤ hZ(h). Therefore τ(h) ≤ Z(h) → 0.

Conversely, suppose τ (h) → 0. For any t ∈ [a, b) and any h ∈ (0, b − t],

x(t + h) − x(t)

= ψ(h, t, x(t)) + τ (h, t).

Taking the limit as h ↓ 0 gives f(t, x(t)) = x

′

(t) = ψ(0, t, x(t)).

Introduction to the Numerical Solution of IVP for ODE

Convergence Theorem for One-Step Methods

Theorem. Suppose f (t, x) is continuous in t, x and uniformly Lipschitz in x on [a, b] × R

Let x(t) be the solution of the IVP x

′

= f (t, x), x(a) = x

on [a, b]. Suppose that the

function ψ(h, t, x) in the one step method satisfies the following two conditions:

1. (Stability) ψ(h, t, x) is continuous in h, t, x and uniformly Lipschitz in x (with Lipschitz

constant K) on 0 ≤ h ≤ h

, a ≤ t ≤ b, x ∈ R

for some h

> 0 with h

≤ b − a, and

2. (Consistency) ψ(0, t, x) = f (t, x).

Let e

(h) = x(t

(h)) − x

(h), where x

is obtained from the one-step method x

i+1

= x

hψ(h, t

, x

). (Note that e

(h) = x

− x

(h) is the error in the initial value x

(h).) Then

(h)| ≤ e

K(t

(h)−a)

(h)| + τ(h)

K(t

(h)−a)

− 1

(h)| ≤ e

K(b−a)

(h)| +

K(b−a)

− 1

τ (h).

Moreover, τ (h) → 0 as h → 0. Therefore, if e

(h) → 0 as h → 0, then

max

0≤i≤

b−a

(h)| → 0 as h → 0,

that is, the approximations converge uniformly on the grid to the solution.

Proof. Hold h > 0 fixed, and ignore rounding error. Subtracting

i+1

= x

+ hψ(h, t

, x

)

from

x(t

i+1

) = x(t

) + hψ(h, t

, x(t

)) + hτ

gives

i+1

| ≤ |e

| + h|ψ(h, t

, x(t

)) − ψ(h, t

, x

)| + h|τ

≤ |e

| + hK|e

| + hτ(h)

= (1 + hK)|e

| + hτ(h).

| ≤ (1 + hK)|e

| + hτ(h),

and

| ≤ (1 + hK)|e

| + hτ(h)

≤ (1 + hK)

| + hτ(h)(1 + (1 + hK)).

By induction,

| ≤ (1 + hK)

| + hτ(h)(1 + (1 + hK) + (1 + hK)

+ · · · + (1 + hK)

i−1

)

= (1 + hK)

| + hτ(h)

(1 + hK)

− 1

(1 + hK) − 1

= (1 + hK)

| + τ(h)

(1 + hK)

− 1

Ordinary Differential Equations

Since (1 + hK)

1
h

↑ e

as h → 0

(for K > 0), and i =

−a

, we have

(1 + hK)

= (1 + hK)

ti−a

≤ e

K(t

−a)

Thus

| ≤ e

K(t

−a)

| + τ(h)

K(t

−a)

− 1

The preceding proposition shows τ (h) → 0, and the theorem follows.

If f is sufficiently smooth, then we know that x(t) ∈ C

p+1

. The theorem thus implies

that if a one-step method is accurate of order p and stable [i.e. ψ is Lipschitz in x], then for
sufficiently smooth f ,

l(h, t) = O(h

p+1

) and thus τ (h) = O(h

If, in addition, e

(h) = O(h

), then

max

(h)| = O(h

i.e. we have p

order convergence of the numerical approximations to the solution.

Example. The “Taylor method of order p” is accurate of order p. If f ∈ C

, then x ∈ C

p+1

and

l(h, t) = x(t + h) −

x(t) + hx

′

(t) + · · · +

(p)

(t)

t+h

(t + h − s)

(p+1)

(s)ds.

|l(h, t)| ≤ M

p+1

(p + 1)!

where M

p+1

= max

a≤t≤b

(p+1)

(t)|.

Fact. A one-step method x

i+1

= x

+ hψ(h, t

, x

) is accurate of order p if and only if

ψ(h, t, x) = T

(h, t, x) + O(h

where T

is the “ψ” for the Taylor method of order p.

Proof. Since

x(t + h) − x(t) = hT

(h, t, x(t)) + O(h

p+1

we have for any given one-step method that

l(h, t) = x(t + h) − x(t) − hψ(h, t, x(t))

= hT

(h, t, x(t)) + O(h

p+1

) − hψ(h, t, x(t))

= h(T

(h, t, x(t)) − ψ(h, t, x(t))) + O(h

p+1

So l(h, t) = O(h

p+1

) iff h(T

− ψ) = O(h

p+1

) iff ψ = T

+ O(h

Remark.

The controlled growth of the effect of the local truncation error (LTE) from previous

steps in the proof of the convergence theorem (a consequence of the Lipschitz continuity of
ψ in x) is called stability. The theorem states:

Stability

Consistency (minimal accuracy)

⇒

Convergence.

In fact, here, the converse is also true.

Introduction to the Numerical Solution of IVP for ODE

Explicit Runge-Kutta methods

One of the problems with Taylor methods is the need to evaluate higher derivatives of f .
Runge-Kutta (RK) methods replace this with the much more reasonable need to evaluate f
more than once to go from x

to x

i+1

. An m-stage (explicit) RK method is of the form

i+1

= x

+ hψ(h, t

, x

with

ψ(h, t, x) =

j=1

(h, t, x),

where a

, . . . , a

are given constants,

(h, t, x) = f (t, x)

and for 2 ≤ j ≤ m,

(h, t, x) = f (t + α

h, x + h

j−1

r=1

(h, t, x))

with α

, . . . , α

and β

(1 ≤ r < j ≤ m) given constants. We usually choose 0 < α

≤ 1,

and for accuracy reasons we choose

(*)

j−1

r=1

(2 ≤ j ≤ m).

Example. m = 2

i+1

= x

+ h(a

(h, t

, x

) + a

(h, t

, x

))

where

(h, t

, x

) = f (t

, x

)

(h, t

, x

) = f (t

+ α

h, x

+ hβ

(h, t

, x

)).

For simplicity, write α for α

and β for β

. Expanding in h,

(h, t, x) = f (t + αh, x + hβf (t, x))

= f (t, x) + αhD

f (t, x) + (D

f (t, x))(hβf (t, x)) + O(h

)

= [f + h(αD

f + β(D

f )f )] (t, x) + O(h

ψ(h, t, x) = (a

+ a

)f + h(a

αD

f + a

β(D

f )f ) + O(h

Recalling that

= f +

f + (D

f )f ),

Ordinary Differential Equations

and that the method is accurate of order two if and only if

ψ = T

+ O(h

we obtain the following necessary and sufficient conditions on a two-stage (explicit) RK
method to be accurate of order two:

+ a

= 1,

α =

1
2

and a

β =

1
2

We require α = β as in (*) (we now see why this condition needs to be imposed), whereupon
these conditions become:

+ a

= 1,

α =

1
2

Therefore, there is a one-parameter family (e.g., parameterized by α) of 2

order, two-stage

(m = 2) explicit RK methods.

Examples.

(1) Setting α =

1
2

gives a

= 1, a

= 0, which is the Modified Euler method.

(2) Choosing α = 1 gives a

1
2

, a

1
2

, which is the Improved Euler method, or Heun’s

method.

Remark.

Note that an m-stage explicit RK method requires m function evaluations (i.e.,

evaluations of f ) in each step (x

to x

i+1

Attainable Orders of Accuracy for Explicit RK methods

# of stages (m)

highest order attainable

← Euler’s method

Explicit RK methods are always stable: ψ inherits its Lipschitz continuity from f .

Example.

Modified Euler

. Let L be the Lipschitz constant for f , and suppose h ≤ h

(for some

≤ b − a).

i+1

= x

+ hf t

, x

f (t

, x

)

ψ(t, h, x) = f t +

, x +

f (t, x)

Introduction to the Numerical Solution of IVP for ODE

|ψ(h, t, x) − ψ(h, t, y)| ≤ L

x +

f (t, x)

− y +

f (t, y)

≤ L|x − y| +

L|f(t, x) − f(t, y)|

≤ L|x − y| +

|x − y|

≤ K|x − y|

where K = L +

is thus the Lipschitz constant for ψ.

Example. A popular 4

order four-stage RK method is

i+1

= x

+ 2k

+ k

)

where

= f (t

, x

)

= f t

, x

= f t

, x

= f (t

+ h, x

+ hk

The same argument as above shows this method is stable.

Remark.

RK methods require multiple function evaluations per step (going from x

to x

i+1

One-step methods discard information from previous steps (e.g., x

i−1

is not used to get x

i+1

— except in its influence on x

). We will next study a class of multi-step methods. But first,

we consider linear difference equations.

Linear Difference Equations (Constant Coefficients)

In this discussion, x

will be a (scalar) sequence defined for i ≥ 0. Consider the linear

difference equation (k-step)

(LDE)

i+k

+ α

k−1

i+k−1

+ · · · + α

= b

(i ≥ 0).

If b

≡ 0, the linear difference equation (LDE) is said to be homogeneous, in which case we

will refer to it as (lh). If b

6= 0 for some i ≥ 0, the linear difference equation (LDE) is said

to be inhomogeneous , in which case we refer to it as (li).

Initial Value Problem (IVP): Given x

for i = 0, . . . , k −1, determine x

satisfying (LDE)

for i ≥ 0.

Theorem. There exists a unique solution of (IVP) for (lh) or (li).

Proof. An obvious induction on i. The equation for i = 0 determines x

, etc.

Theorem. The solution set of (lh) is a k-dimensional vector space (a subspace of the set of
all sequences {x

}

i≥0

Ordinary Differential Equations

Proof Sketch. Choosing







...

k−1







= e

∈ R

for j = 1, 2, . . . k and then solving (lh) gives a basis of the solution space of (lh).

Define the characteristic polynomial of (lh) to be

p(r) = r

+ α

k−1

+ · · · + α

Let us assume that α

6= 0. (If α

= 0, (LDE) isn’t really a k-step difference equation since

we can shift indices and treat it as a e

k-step difference equation for a e

k < k, namely e

k = k −ν,

where ν is the smallest index with α

6= 0.) Let r

, . . . , r

be the distinct zeroes of p, with

multiplicities m

, . . . , m

. Note that each r

6= 0 since α

6= 0, and m

+ · · · + m

= k. Then

a basis of solutions of (lh) is:

}

∞

i=0

: 1 ≤ j ≤ s, 0 ≤ l ≤ m

− 1

Example. Fibonacci Sequence:

i+2

− F

i+1

− F

= 0,

= 1.

The associated characteristic polynomial r

− r − 1 = 0 has roots

1 ±

√

≈ 1.6, r

−

≈ −0.6).

The general solution of (lh) is

= C

1 +

√

+ C

−

1 −

√

The initial conditions F

= 0 and F

= 1 imply that C

√

and C

−

= −

√

. Hence

√





1 +

√

−

1 −

√



 .

Since |r

−

| < 1, we have

1 −

√

→ 0 as i → ∞.

Hence, the Fibonacci sequence behaves asymptotically like the sequence

√

Introduction to the Numerical Solution of IVP for ODE

Remark.

If α

= α

= · · · = α

ν−1

= 0 and α

6= 0 (i.e., 0 is a root of multiplicity ν), then

, x

, . . . , x

ν−1

are completely independent of x

for i ≥ ν. So x

i+k

+ · · · + α

i+ν

= b

for

i ≥ 0 with x

given for i ≥ ν behaves like a (k − ν)-step difference equation.

Remark.

Define e







i+1

i+k−1







. Then e

i+1

= Ae

for i ≥ 0, where

A =







. .. ...

−α

. . .

−α

k−1







and e







...

k−1







is given by the I.C. So (lh) is equivalent to the one-step vector difference

equation

i+1

= Ae

i ≥ 0,

whose solution is e

= A

. The characteristic polynomial of (lh) is the characteristic

polynomial of A. Because A is a companion matrix, each distinct eigenvalue has only one
Jordan block. If A = P JP

−1

is the Jordan decomposition of A (J in Jordan form, P

invertible), then

= P J

−1

Let J

be the m

× m

block corresponding to r

(for 1 ≤ j ≤ s), so J

= r

I + Z

, where Z

denotes the m

× m

shift matrix:







. .. ...

1
0







Then

= (r

I + Z

)

l=0

i
l

i−l

Since

i
l

is a polynomial in i of degree l and Z

= 0, we see entries in e

of the form

(constant) i

for 0 ≤ l ≤ m

− 1.

Remark.

(li) becomes

i+1

= Ae

+ eb

i ≥ 0,

Ordinary Differential Equations

where eb

= [0, . . . , 0, b

]

. This leads to a discrete version of Duhamel’s principle (exercise).

Remark.

All solutions {x

}

i≥0

of (lh) stay bounded (i.e. are elements of l

∞

)

⇔ the matrix A is power bounded (i.e., ∃ M so that kA

k ≤ M for all i ≥ 0)

⇔ the Jordan blocks J

, . . . , J

are all power bounded

⇔

(a) each |r

| ≤ 1

(for 1 ≤ j ≤ s)

and (b) for any j with m

> 1 (multiple roots), |r

| < 1

If (a) and (b) are satisfied, then the map e

7→ {x

}

i≥0

is a bounded linear operator from R

(or C

) into l

∞

(exercise).

Linear Multistep Methods (LMM)

A LMM is a method of the form

j=0

i+j

= h

j=0

i+j

i ≥ 0

for the approximate solution of an ODE IVP

′

= f (t, x),

x(a) = x

Here we want to approximate the solution x(t) of this IVP for a ≤ t ≤ b at the points
t

= a + ih (where h is the time step), 0 ≤ i ≤

b−a

. The term x

denotes the approximation

to x(t

). We have set f

i+j

= f (t

i+j

, x

i+j

). We normalize the coefficients so that α

= 1.

The above is called a k-step LMM (if at least one of the coefficients α

and β

is non-zero).

The above equation is similar to a difference equation in that one is solving for x

i+k

from

, x

i+1

, . . . , x

i+k−1

. We assume as usual that f is continuous in (t, x) and uniformly Lipschitz

in x. For simplicity of notation, we will assume that x(t) is real and scalar; the analysis that
follows applies to x(t) ∈ R

or x(t) ∈ C

(viewed as R

for differentiability) with minor

adjustments.

Example. (Midpoint rule) (explicit)

x(t

i+2

) − x(t

) =

i+2

′

(s)ds ≈ 2hx

′

i+1

) = 2hf (t

i+1

, x(t

i+1

)).

This approximate relationship suggests the LMM

Midpoint rule

i+2

− x

= 2hf

i+1

Example. (Trapezoid rule) (implicit)

The approximation

x(t

i+1

) − x(t

) =

i+1

′

(s)ds ≈

′

i+1

) + x

′

))

Introduction to the Numerical Solution of IVP for ODE

suggests the LMM

Trapezoid rule

i+1

− x

i+1

+ f

) .

Explicit vs Implicit.

If β

= 0, the LMM is called explicit: once we know x

, x

i+1

, . . . , x

i+k−1

, we compute imme-

diately

i+k

k−1

j=0

(hβ

i+j

− α

i+j

) .

On the other hand, if β

6= 0, the LMM is called implicit: knowing x

, x

i+1

, . . . , x

i+k−1

, we

need to solve

i+k

= hβ

f (t

i+k

, x

i+k

) −

k−1

j=0

(α

i+j

− hβ

i+j

)

for x

i+k

Remark.

If h is sufficiently small, implicit LMM methods also have unique solutions given h

and x

, x

, . . . , x

k−1

. To see this, let L be the Lipschitz constant for f . Given x

, . . . , x

i+k−1

the value for x

i+k

is obtained by solving the equation

i+k

= hβ

f (t

i+k

, x

i+k

) + g

where

k−1

j=0

(hβ

i+j

− α

i+j

)

is a constant as far as x

i+k

is concerned. That is, we are looking for a fixed point of

Φ(x) = hβ

f (t

i+k

, x) + g

Note that if h|β

|L < 1, then Φ is a contraction:

|Φ(x) − Φ(y)| ≤ h|β

| |f(t

i+k

, x) − f(t

i+k

, y)| ≤ h|β

|L|x − y|.

So by the Contraction Mapping Fixed Point Theorem, Φ has a unique fixed point. Any initial
guess for x

i+k

leads to a sequence converging to the fixed point using functional iteration

(l+1)
i+k

= hβ

f (t

i+k

, x

(l)
i+k

) + g

which is initiated at some initial point x

(0)
i+k

. In practice, one chooses either

(1) iterate to convergence, or

(2) a fixed number of iterations, using an explicit method to get the initial guess x

(0)
i+k

This pairing is often called a Predictor-Corrector Method.

Ordinary Differential Equations

Function Evaluations. One FE means evaluating f once.

Explicit LMM: 1 FE per step (after initial start)
Implicit LMM: ? FEs per step if iterate to convergence

usually 2 FE per step for a Predictor-Corrector Method.

Initial Values. To start a k-step LMM, we need x

, x

, . . . , x

k−1

. We usually take x

= x

and approximate x

, . . . , x

k−1

using a one-step method (e.g., a Runge-Kutta method).

Local Truncation Error. For a true solution x(t) to x

′

= f (t, x), define the LTE to be

l(h, t) =

j=0

x(t + jh) − h

j=0

′

(t + jh).

If x ∈ C

p+1

, then

x(t + jh) = x(t) + jhx

′

(t) + · · · +

(jh)

(p)

(t) + O(h

p+1

) and

′

(t + jh) = hx

′

(t) + jh

′′

(t) + · · · +

p−1

(p − 1)!

(p)

(t) + O(h

p+1

)

and so

l(h, t) = C

x(t) + C

′

(t) + · · · + C

(p)

(t) + O(h

p+1

where

= α

+ · · · + α

= (α

+ 2α

+ · · · + kα

) − (β

+ · · · + β

)

...

(α

+ 2

+ · · · + k

) −

(q − 1)!

(β

+ 2

q−1

+ · · · + k

q−1

Definition. A LMM is called accurate of order p if l(h, t) = O(h

p+1

) for any solution of

′

= f (t, x) which is C

p+1

Fact. A LMM is accurate of order at least p iff C

= C

= · · · = C

= 0. (It is called

accurate of order exactly p if also C

p+1

6= 0.)

Remarks.

(i) For the LTE of a method to be o(h) for all f ’s, we must have C

= C

= 0. To see

this, for any f which is C

, all solutions x(t) are C

, so

l(h, t) = C

x(t) + C

′

(t) + O(h

) is o(h) iff C

= C

= 0 .

(ii) Note that C

, C

, . . . depend only on α

, . . . , α

, β

, . . . , β

and not on f . So here,

“minimal accuracy” is first order.

Introduction to the Numerical Solution of IVP for ODE

Definition. A LMM is called consistent if C

= C

= 0 (i.e., at least first-order accurate).

Remark.

If a LMM is consistent, then any solution x(t) for any f (continuous in (t, x),

Lipschitz in x) has l(h, t) = o(h). To see this, note that since x ∈ C

x(t + jh) = x(t) + jhx

′

(t) + o(h) and hx

′

(t + jh) = hx

′

(t) + o(h),

l(h, t) = C

x(t) + C

′

(t) + o(h).

Exercise: Verify that the o(h) terms converge to 0 uniformly in t (after dividing by h) as
h → 0: use the uniform continuity of x

′

(t) on [a, b].

Definition. A k-step LMM

i+j

= h

i+j

is called convergent if for each IVP x

′

= f (t, x), x(a) = x

on [a, b] (f ∈ (C, Lip)) and for

any choice of x

(h), . . . , x

k−1

(h) for which

lim

h→0

|x(t

(h)) − x

(h)| = 0 for i = 0, . . . , k − 1,

we have

lim

h→0

max

{i:a≤t

(h)≤b}

|x(t

(h)) − x

(h)| = 0 .

Remarks.

(i) This asks for uniform decrease of the error on the grid as h → 0.

(ii) By continuity of x(t), the condition on the initial values is equivalent to x

(h) → x

for i = 0, 1, . . . , k − 1.

Fact. If a LMM is convergent, then the zeroes of the (first) characteristic polynomial of the
method p(r) = α

+ · · · + α

satisfy the Dahlquist root condition:

(a) all zeroes r satisfy |r| ≤ 1, and

(b) multiple zeroes r satisfy |r| < 1.

Examples. Consider the IVP x

′

= 0, a ≤ t ≤ b, x(a) = 0. So f ≡ 0. Consider the LMM:

i+j

= 0 .

Ordinary Differential Equations

(1) Let r be any zero of p(r). Then the solution with initial conditions

= hr

for 0 ≤ i ≤ k − 1

= hr

for 0 ≤ i ≤

b − a

Suppose h =

b−a

for some m ∈ Z. If the LMM is convergent, then

(h) → x(b) = 0

as m → ∞. But

(h) = hr

b − a

(h) − x(b)| =

b − a

| → 0 as m → ∞

iff |r| ≤ 1.

(2) Similarly if r is a multiple zero of p(r), taking x

(h) = hir

for 0 ≤ i ≤ k − 1 gives

(h) = hir

0 ≤ i ≤

b − a

So if h =

b−a

, then

(h) =

b − a

= (b − a)r

so x

(h) → 0 as h → 0 iff |r| < 1.

Definition. A LMM is called zero-stable if it satisfies the Dahlquist root condition.

Recall from our discussion of linear difference equations that zero-stability is equivalent to
requiring that all solutions of (lh)

k
j=0

i+j

= 0 for i ≥ 0 stay bounded as i → ∞.

Remark.

A consistent one-step LMM (i.e., k = 1) is always zero-stable. Indeed, consistency

implies that C

= C

= 0, which in turn implies that p(1) = α

+α

= C

= 0 and so r = 1 is

the

zero of p(r). Therefore α

= 1, α

= −1, so the characteristic polynomial is p(r) = r − 1,

and the LMM is zero-stable.

Exercise: Show that if an LMM is convergent, then it is consistent.

Key Theorem. [LMM Convergence]

A LMM is convergent if and only if it is zero-stable and consistent. Moreover, for zero-stable
methods, we get an error estimate of the form

max

a≤t

(h)≤b

|x(t

(h)) − x

(h)| ≤ K

max

0≤i≤k−1

|x(t

(h)) − x

(h)|

}

initial error

max

|l(h, t

(h))|

}

“growth of error”

controlled by

zero-stability

Introduction to the Numerical Solution of IVP for ODE

Remark.

If x ∈ C

p+1

and the LMM is accurate of order p, then |LT E|/h = O(h

). To get

-order convergence (i.e., LHS = O(h

)), we need

(h) = x(t

(h)) + O(h

) for i = 0, . . . , k − 1.

This can be done using k − 1 steps of a RK method of order ≥ p − 1.

Lemma. Consider

(li)

j=0

i+j

= b

for i ≥ 0

(where α

= 1),

with the initial values x

, . . . , x

k−1

given, and suppose that the characteristic polynomial

p(r) =

k
j=0

satisfies the Dahlquist root condition. Then there is an M > 0 such that

for i ≥ 0,

i+k

| ≤ M

max{|x

|, . . . , |x

k−1

|} +

ν=0

Remark.

Recall that the Dahlquist root condition implies that there is an M > 0 for which

∞

≤ M for all i ≥ 0, where

A =







. .. ...

−α

· · ·

−α

k−1







is the companion matrix for p(r), and k · k

∞

is the operator norm induced by the vector

norm k · k

∞

. The M in the Lemma can be taken to be the same as this M bounding kA

∞

Proof. Let e

= [x

, x

i+1

, . . . , x

i+k−1

]

and eb

= [0, . . . , 0, b

]

. Then e

i+1

= Ae

+ eb

, so

by induction

i+1

= A

i+1

ν=0

i−ν

Thus

i+k

| ≤ ke

i+1

∞

≤ kA

i+1

∞

ν=0

i−ν

∞

keb

∞

≤ M(ke

∞

ν=0

|).

Ordinary Differential Equations

Proof of the LMM Convergence Theorem. The fact that convergence implies zero-
stability and consistency has already been discussed. Suppose a LMM is zero-stable and
consistent. Let x(t) be the true solution of the IVP x

′

= f (t, x), x(a) = x

on [a, b], let L be

the Lipschitz constant for f , and set

β =

j=0

|β

Hold h fixed, and set

(h) = x(t

(h)) − x

(h),

E = max{|e

|, . . . , |e

k−1

|},

(h) = l(h, t

(h)),

λ(h) = max

i∈I

(h)|,

where I = {i ≥ 0 : i + k ≤

b−a

Step 1. The first step is to derive a “difference inequality” for |e

|. This difference inequality

is a discrete form of the integral inequality leading to Gronwall’s inequality. For i ∈ I, we
have

j=0

x(t

i+j

) = h

j=0

f (t

i+j

, x(t

i+j

)) + l

j=0

i+j

= h

j=0

i+j

Subtraction gives

j=0

i+j

= b

where

≡ h

j=0

(f (t

i+j

, x(t

i+j

)) − f(t

i+j

, x

i+j

)) + l

Then

| ≤ h

j=0

|β

|L|e

i+j

| + |l

So, by the preceeding Lemma with x

i+k

replaced by e

i+k

, we obtain for i ∈ I

i+k

| ≤ M

E +

ν=0

≤ M

E + hL

ν=0

j=0

|β

||e

ν+j

| +

ν=0

≤ M

E + hL|β

||e

i+k

| + hLβ

i+k−1

ν=0

| +

ν=0

Introduction to the Numerical Solution of IVP for ODE

From here on, assume h is small enough that

MhL|β

| ≤

1
2

(Since

h ≤ b − a : MhL|β

| ≥

1
2

is a compact subset of (0, b − a], the estimate in the Key

Theorem is clearly true for those values of h.) Moving MhL|β

||e

i+k

| to the LHS gives

i+k

| ≤ hM

i+k−1

ν=0

| + M

E + M

λ/h

for i ∈ I, where M

= 2MLβ, M

= 2M, and M

= 2M(b−a). (Note: For explicit methods,

= 0, so the restriction MhL|β

| ≤

1
2

is unnecessary, and the factors of 2 in M

, M

can be dropped.)

Step 2. We now use a discrete “comparison” argument to bound |e

|. Let y

be the solution

(∗)

i+k

= hM

i+k−1

ν=0

+ (M

E + M

λ/h) for i ∈ I,

with initial values y

= |e

| for 0 ≤ j ≤ k − 1. Then clearly by induction |e

i+k

| ≤ y

i+k

for

i ∈ I. Now

≤ hM

kE + (M

E + M

λ/h) ≤ M

E + M

λ/h,

where M

= (b − a)M

k + M

. Subtracting (∗) for i from (∗) for i + 1 gives

i+k+1

− y

i+k

= hM

i+k

and so y

i+k+1

= (1 + hM

i+k

Therefore, by induction we obtain for i ∈ I:

i+k

= (1 + hM

)

≤ (1 + hM

)

(b−a)/h

≤ e

(b−a)

≤ K

E + K

λ/h,

where K

= e

(b−a)

and K

= e

(b−a)

. Thus, for i ∈ I,

i+k

| ≤ K

E + K

λ/h;

since K

≥ M

≥ M ≥ 1, also |e

| ≤ E ≤ K

E + K

λ/h for 0 ≤ j ≤ k − 1. Since

consistency implies λ = o(h), we are done.

Remarks.

(1) Note that K

and K

depend only on a, b, L, k, the α

’s and β

’s, and M.

(2) The estimate can be refined — we did not try to get the best constants K

, K

. For

example, e

(b−a)

could be replaced by e

−a)

in both K

and K

, yielding more

precise estimates depending on i, similar to the estimate for one-step methods.

Wyszukiwarka

Podobne podstrony:
Bradley Numerical Solutions of Differential Equations [sharethefiles com]
Accelerating numerical solution of Stochastic DE with CUDA
Numerical Analysis of Conditions for Ignition of Compact Metal Specimens and Foil in Oxygen
An approximate solution of the Riemann problem for a realisable second moment turbulent closure
Barret et al Templates For The Solution Of Linear Systems Building Blocks For Iterative Methods [s
37 509 524 Microstructure and Wear Resistance of HSS for Rolling Mill Rolls
Preparation of Material for a Roleplaying?venture
fema declaration of lack of workload for pr npsc 2008
Overview of ODEs for Computational Science
76 1075 1088 The Effect of a Nitride Layer on the Texturability of Steels for Plastic Moulds
45 625 642 Numerical Simulation of Gas Quenching of Tool Steels
Coupling of Technologies for Concurrent ECD and Barite Sag Management
37 509 524 Microstructure and Wear Resistance of HSS for Rolling Mill Rolls
Majewski, Marek; Bors, Dorota On the existence of an optimal solution of the Mayer problem governed
The Fiqh of Hajj for Women
variation of 7979 for claire opis
Wicca Book of Spells and Witchcraft for Beginners The Guide of Shadows for Wiccans, Solitary Witche
The Treasure of Treasures for Alchemists by Paracelsus
Dynkin E B Superdiffusions and Positive Solutions of Nonlinear PDEs (web draft,2005)(108s) MCde

więcej podobnych podstron