IVP of ODE

(1)

IVP of ODE

¹

Initial-Value Problems for Ordinary Differential Equations

NTNU

Tsung-Min Hwang December 8, 2003

Department of Mathematics – NTNU Tsung-Min Hwang December 8, 2003

(2)

IVP of ODE

²

1 – Existence and Uniqueness of Solutions . . . . 4

2 – Euler’s Method . . . . 6

3 – Runge-Kutta Methods . . . . 9

4 – Multistep Methods . . . 15

5 – Systems and Higher-Order Ordinary Differential Equations . . . 25

(3)

IVP of ODE

³

In this chapter, we discuss numerical methods for solving ordinary differential equations of initial-value problems (IVP) of the form







y⁰ = f (x, y), x ∈ [a, b]

y(x₀) = y₀, ⁽¹⁾

where y is a function of x^, f is a function of y ^and x^, x₀ is called the initial point, and y₀ ^the

initial value. The numerical values of y(x) on an interval containing x₀ ^{are to be}

determined.

(4)

IVP of ODE

⁴

1 – Existence and Uniqueness of Solutions

Theorem 1 If f (x, y) is continuous in a region Ω^{, where}

Ω = {(x, y); |x − x₀| ≤ α, |y − y₀| ≤ β} ⁽²⁾

then the IVP (1) has a solution y(x) ^for |x − x₀| ≤ min{α, β

M }^{, where} M = max

(x,y)∈Ω |f (x, y)|^.

Theorem 2 If f ^and ^∂f_∂x are continuous in Ω, then the IVP (1) has a unique solution in the interval |x − x₀| ≤ min{α, _M^β }^.

Theorem 3 If f is continuous in a ≤ x ≤ b^, −∞ < y < ∞ ^and

|f (x, y₁) − f (x, y₂)| ≤ L|y₁ − y₂|

for some positive constant L, (that is, f is Lipschitz continuous in y), then IVP (1) has a unique solution in the interval [a, b]^.

(5)

IVP of ODE

⁴

Ω = {(x, y); |x − x₀| ≤ α, |y − y₀| ≤ β} ⁽²⁾

(x,y)∈Ω |f (x, y)|^.

|f (x, y₁) − f (x, y₂)| ≤ L|y₁ − y₂|

(6)

IVP of ODE

⁴

Ω = {(x, y); |x − x₀| ≤ α, |y − y₀| ≤ β} ⁽²⁾

(x,y)∈Ω |f (x, y)|^.

|f (x, y₁) − f (x, y₂)| ≤ L|y₁ − y₂|

(7)

IVP of ODE

⁵

☞ Numerical integration of ODEs : Given

dy

dx = f (x, y), y(x₀) = y₀.

Find approximate solution at discrete values of x xj = x₀ + jh

where h is the “stepsize”.

☞ Graphical interpretation

(8)

IVP of ODE

⁶

2 – Euler’s Method

On of the simplest methods for solving the IVP (1) is the Euler’s method.

Subdivide [a, b] ^into n subintervals of equal length h with mesh points {x₀, x₁, . . . , x_n}

where

xi = a + ih, ∀ i = 0, 1, 2, . . . , n.

Recall the Taylor’s Theorem

y(x_i+1) = y(xi) + (x_i+1 − xi)y⁰(xi) + (x_i+1 − xi)²

2 y⁰⁰(ξi)

= y(x_i) + hy⁰(x_i) + h²

2 y⁰⁰(ξ_i) ⁽³⁾

= y(x_i) + hf (x_i, y(x_i)) + h²

2 y⁰⁰(ξ_i)

for some ξi ∈ [xi, xi+1]^.

(9)

IVP of ODE

⁶

where

xi = a + ih, ∀ i = 0, 1, 2, . . . , n.

2 y⁰⁰(ξi)

= y(x_i) + hy⁰(x_i) + h²

2 y⁰⁰(ξ_i) ⁽³⁾

= y(x_i) + hf (x_i, y(x_i)) + h²

2 y⁰⁰(ξ_i)

(10)

IVP of ODE

⁶

where

xi = a + ih, ∀ i = 0, 1, 2, . . . , n.

2 y⁰⁰(ξi)

= y(x_i) + hy⁰(x_i) + h²

2 y⁰⁰(ξ_i) ⁽³⁾

= y(x_i) + hf (x_i, y(x_i)) + h²

2 y⁰⁰(ξ_i)

(11)

IVP of ODE

⁷

We have the formulation of Euler’s mehtod

x_k+1 = xk + h = x₀ + (k + 1)h y_k+1 = yk + hf (xk, yk)

On the other hand, from (3),

f (xi, y(xi)) = y⁰(xi) = y(x_i+1) − y(xi)

h + h

2y⁰⁰(ξi)

It implies that Euler’s method is a first-order (O(h)) ^method.

The Improved Euler’s method

x_k+1 = xk + h = x₀ + (k + 1)h

y_k^∗₊₁ = yk + hf (xk, yk) ⁽⁴⁾

yk+1 = yk + h

2 f(xk, yk) + f (xk+1, y_k^∗₊₁)

(12)

IVP of ODE

⁷

h + h

2y⁰⁰(ξi)

x_k+1 = xk + h = x₀ + (k + 1)h

y_k^∗₊₁ = yk + hf (xk, yk) ⁽⁴⁾

yk+1 = yk + h

2 f(xk, yk) + f (xk+1, y_k^∗₊₁)

(13)

IVP of ODE

⁷

h + h

2y⁰⁰(ξi)

x_k+1 = xk + h = x₀ + (k + 1)h

y_k^∗₊₁ = yk + hf (xk, yk) ⁽⁴⁾

yk+1 = yk + h

2 f(xk, yk) + f (xk+1, y_k^∗₊₁)

(14)

IVP of ODE

⁷

h + h

2y⁰⁰(ξi)

x_k+1 = xk + h = x₀ + (k + 1)h

y_k^∗₊₁ = yk + hf (xk, yk) ⁽⁴⁾

yk+1 = yk + h

2 f(xk, yk) + f (xk+1, y_k^∗₊₁)

(15)

IVP of ODE

⁸

Example 1 Use Euler’s method to integrate

dy

dx = x − 2y, y(0) = 1.

The exact solution is

y = 1

4 2x − 1 + 5e⁻^2x .

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

y

Example 1 for Euler method

Exact h = 0.2 h = 0.1 h = 0.05

(16)

IVP of ODE

⁸

dy

dx = x − 2y, y(0) = 1.

y = 1

4 2x − 1 + 5e⁻^2x .

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

y

Exact h = 0.2 h = 0.1 h = 0.05

(17)

IVP of ODE

⁸

dy

dx = x − 2y, y(0) = 1.

y = 1

4 2x − 1 + 5e⁻^2x .

0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x

y

Exact h = 0.2 h = 0.1 h = 0.05

(18)

IVP of ODE

⁹

3 – Runge-Kutta Methods

One of the most important methods for solving the IVP (1) is the Runge-Kutta method.

y(x + h) = y(x) + hy⁰(x) + h²

2 y⁰⁰ + h³

2 y⁰⁰⁰ + . . .

= y(x) + hy⁰(x) + h²

2 y⁰⁰ + O(h³)

(19)

IVP of ODE

⁹

y(x + h) = y(x) + hy⁰(x) + h²

2 y⁰⁰ + h³

2 y⁰⁰⁰ + . . .

= y(x) + hy⁰(x) + h²

2 y⁰⁰ + O(h³)

(20)

IVP of ODE

⁹

y(x + h) = y(x) + hy⁰(x) + h²

2 y⁰⁰ + h³

2 y⁰⁰⁰ + . . .

= y(x) + hy⁰(x) + h²

2 y⁰⁰ + O(h³)

(21)

IVP of ODE

¹⁰

By differentiating y(x)^{, we have}

y⁰ = f (x, y) ≡ f

y⁰⁰ = d

dxf = fx + fyy⁰ = fx + fyf y⁰⁰⁰ = fxx + fxyf + f d

dxfy + fy

d dxf

= fxx + fxyf + f (f_y dy

dy

dx + fyx) + fy(fx + fyf )

= fxx + fxyf + f (fyx + fyyf ) + fy(fx + fyf )

(22)

IVP of ODE

¹⁰

y⁰ = f (x, y) ≡ f y⁰⁰ = d

dxf = fx + fyy⁰ = fx + fyf

y⁰⁰⁰ = fxx + fxyf + f d

dxfy + fy

d dxf

dy

(23)

IVP of ODE

¹⁰

y⁰ = f (x, y) ≡ f y⁰⁰ = d

dxfy + fy

d dxf

dy

(24)

IVP of ODE

¹⁰

y⁰ = f (x, y) ≡ f y⁰⁰ = d

dxfy + fy

d dxf

dy

(25)

IVP of ODE

¹⁰

y⁰ = f (x, y) ≡ f y⁰⁰ = d

dxfy + fy

d dxf

dy

(26)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

Apply Taylor’s Theorem on f

f (x + h, y + hf ) = f (x, y) + hfx + hf fy + O(h²)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(27)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(28)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(29)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(30)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(31)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(32)

IVP of ODE

¹¹

Then

y(x + h) = y + hy⁰ + h²

2 y⁰⁰ + O(h³)

= y + hf + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2f + h²

2 (fx + fyf ) + O(h³)

= y + h

2f + h

2(f + hfx + hfyf ) + O(h³)

⇒ f + hfx + hf fy ≈ f (x + h, y + hf )

⇒ y(x + h) = y + 1

2hf + h

2f (x + h, y + hf ) + O(h³)

(33)

IVP of ODE

¹²

Let

F₁ = hf (x, y), F₂ = hf (x + h, y + F₁)

Then

y(x + h) = y + 1

2(F₁ + F₂).

This is called the second-order Runge-Kutta method. Two function evaluations are required at each step.

Algorithm 1 Algorithm for second-order Runge-Kutta mehtod : For k = 0, 1, 2, . . .

xk+1 = xk + h = x₀ + (k + 1)h F₁ = hf (xk, yk)

F₂ = hf (x_k+1, y_k + F₁) yk+1 = yk + ¹₂(F₁ + F₂)

End for Department of Mathematics – NTNU Tsung-Min Hwang December 8, 2003

(34)

IVP of ODE

¹²

Let

F₁ = hf (x, y), F₂ = hf (x + h, y + F₁)

Then

y(x + h) = y + 1

2(F₁ + F₂).

F₂ = hf (x_k+1, y_k + F₁) yk+1 = yk + ¹₂(F₁ + F₂)

(35)

IVP of ODE

¹²

Let

F₁ = hf (x, y), F₂ = hf (x + h, y + F₁)

Then

y(x + h) = y + 1

2(F₁ + F₂).

This is called the second-order Runge-Kutta method.

Two function evaluations are required at each step.

F₂ = hf (x_k+1, y_k + F₁) yk+1 = yk + ¹₂(F₁ + F₂)

(36)

IVP of ODE

¹²

Let

F₁ = hf (x, y), F₂ = hf (x + h, y + F₁)

Then

y(x + h) = y + 1

2(F₁ + F₂).

F₂ = hf (x_k+1, y_k + F₁) yk+1 = yk + ¹₂(F₁ + F₂)

(37)

IVP of ODE

¹²

Let

F₁ = hf (x, y), F₂ = hf (x + h, y + F₁)

Then

y(x + h) = y + 1

2(F₁ + F₂).

F₂ = hf (x_k+1, y_k + F₁) yk+1 = yk + ¹₂(F₁ + F₂)

End for

(38)

IVP of ODE

¹³

General form of second-order Runge-Kutta method :

y(x + h) = y + ω₁hf (x, y) + ω₂hf (x + αh, y + βhf ) + O(h³)

where ω₁, ω₂, α, β are constants to be defined, and

ω₁ + ω₂ = 1, ω₂α = 1

2, ω₂β = 1 2

By letting ω₁ = 0, ω₂ = 1, α = β = ¹₂ leads to the modified Euler’s mehtod. Fourth-Order Runge-kutta method

y(x + h) = y + 1

6(F₁ + 2F₂ + 2F₃ + F₄) + O(h⁵)

where











F₁ = hf (x, y)

F₂ = hf (x + ¹₂h, y + ¹₂F₁) F₃ = hf (x + ¹₂h, y + ¹₂F₂) F₄ = hf (x + h, y + F₃)

(39)

IVP of ODE

¹³

ω₁ + ω₂ = 1, ω₂α = 1

2, ω₂β = 1 2

By letting ω₁ = 0, ω₂ = 1, α = β = ¹₂ leads to the modified Euler’s mehtod.

Fourth-Order Runge-kutta method

y(x + h) = y + 1

6(F₁ + 2F₂ + 2F₃ + F₄) + O(h⁵)

where











F₁ = hf (x, y)

(40)

IVP of ODE

¹³

ω₁ + ω₂ = 1, ω₂α = 1

2, ω₂β = 1 2

By letting ω₁ = 0, ω₂ = 1, α = β = ¹₂ leads to the modified Euler’s mehtod.

Fourth-Order Runge-kutta method

y(x + h) = y + 1

6(F₁ + 2F₂ + 2F₃ + F₄) + O(h⁵)

where











F₁ = hf (x, y)

(41)

IVP of ODE

¹⁴

Algorithm 2 Algorithm for fourth-order Runge-Kutta mehtod : For k = 0, 1, 2, . . .

x_k+¹

2 = xk + ¹₂h

x_k+1 = xk + h = x₀ + (k + 1)h F₁ = hf (x_k, y_k)

F₂ = hf (x_k+¹

2, yk + ¹₂F₁) F₃ = hf (x_k+¹

2, yk + ¹₂F₂) F₄ = hf (x_k+1, y_k + F₃)

y_k+1 = yk + ¹₆(F₁ + 2F₂ + 2F₃ + F₄)

End for

The fourth Runge-Kutta method involves a local truncation error of O(h⁵)^.

(42)

IVP of ODE

¹⁵

4 – Multistep Methods

Example 2 (Stiff problem, van der Pol equation)

y₁⁰ = y₂,

y₂⁰ = µ(1 − y₁²)y₂ − y₁.

0 500 1000 1500 2000 2500 3000

−2.5

−2

−1.5

−1

−0.5 0 0.5 1 1.5 2 2.5

Solution of van der Pol Equation, µ = 1000

time t solution y 1

(43)

IVP of ODE

¹⁶

805.715 805.72 805.725 805.73 805.735

−2

−1.5

−1

−0.5 0 0.5

Solution of van der Pol Equation, µ = 1000

time t solution y 1

(44)

IVP of ODE

¹⁷

Definition 1 A multistep method for solving the initial-value problem

y⁰ = f (x, y), a ≤ x ≤ b, y(a) = α

is

w_i+1 = a_m−1wi + a_m−2w_i−1 + · · · + a₀w_i+1−m +

h [bmf (x_i+1, w_i+1) + b_m−1f (xi, wi) + · · · + b₀f (x_i+1−m, w_i+1−m)]

for i = m − 1, m, . . . , N − 1, where the starting values

w₀ = α, w₁ = α₁, · · · , wm−1 = αm−1

are specified and h = (b − a)/N^.

☞ If b_m = 0, the method is called explicit or open

☞ _If bm 6= 0, the method is called implicit or closed. Department of Mathematics – NTNU Tsung-Min Hwang December 8, 2003

(45)

IVP of ODE

¹⁷

y⁰ = f (x, y), a ≤ x ≤ b, y(a) = α

is

w_i+1 = a_m−1wi + a_m−2w_i−1 + · · · + a₀w_i+1−m +

w₀ = α, w₁ = α₁, · · · , wm−1 = αm−1

☞ _If bm 6= 0, the method is called implicit or closed. Department of Mathematics – NTNU Tsung-Min Hwang December 8, 2003

(46)

IVP of ODE

¹⁷

y⁰ = f (x, y), a ≤ x ≤ b, y(a) = α

is

w_i+1 = a_m−1wi + a_m−2w_i−1 + · · · + a₀w_i+1−m +

w₀ = α, w₁ = α₁, · · · , wm−1 = αm−1

☞ _If bm 6= 0, the method is called implicit or closed.

(47)

IVP of ODE

¹⁸

Example 3

☞ Explicit fourth-order Adams-Bashforth method:

w₀ = α, w₁ = α₁, w₂ = α₂, w₃ = α₃, w_i+1 = w_i + h

24 [55f (x_i, w_i) − 59f (x_i−1, w_i−1) +37f (xi−2, wi−2) − 9f (xi−3, wi−3)]

for each i = 3, 4, . . . , N − 1^.

☞ Implicit fourth-order Adams-Moulton method:

w₀ = α, w₁ = α₁, w₂ = α₂, w_i+1 = w_i + h

24 [9f (x_i+1, w_i+1) + 19f (x_i, w_i)

−5f (x_i−1, w_i−1) + f (x_i−2, w_i−2)]

for each i = 2, 3, . . . , N − 1^.

(48)

IVP of ODE

¹⁸

Example 3

☞ Explicit fourth-order Adams-Bashforth method:

w₀ = α, w₁ = α₁, w₂ = α₂, w₃ = α₃, w_i+1 = w_i + h

24 [55f (x_i, w_i) − 59f (x_i−1, w_i−1) +37f (xi−2, wi−2) − 9f (xi−3, wi−3)]

for each i = 3, 4, . . . , N − 1^.

☞ Implicit fourth-order Adams-Moulton method:

w₀ = α, w₁ = α₁, w₂ = α₂, w_i+1 = w_i + h

24 [9f (x_i+1, w_i+1) + 19f (x_i, w_i)

−5f (x_i−1, w_i−1) + f (x_i−2, w_i−2)]

for each i = 2, 3, . . . , N − 1^.

(49)

IVP of ODE

¹⁹

➳ Derivation of the multistep methods:

y(xi+1) − y(xi) =

Z x_i+1 xi

y⁰(x)dx =

Z x_i+1 xi

f (x, y(x))dx

⇒ y(x_i+1) = y(xi) +

Z xi+1

xi

f (x, y(x))dx

P (x): an interpolating polynomial to f (x, y(x)) ^and y(x_i) ≈ w_i^{. Then} y(x_i+1) ≈ w_i +

Z xi+1

xi

P (x)dx

➳ Adams-Bashforth explicit m-step technique: The binomial formula is defined as

s k

= s(s − 1) · · · (s − k + 1) k!

(50)

IVP of ODE

¹⁹

y(xi+1) − y(xi) =

Z x_i+1 xi

y⁰(x)dx =

Z x_i+1 xi

f (x, y(x))dx

⇒ y(x_i+1) = y(xi) +

Z xi+1

xi

f (x, y(x))dx

P (x): an interpolating polynomial to f (x, y(x)) ^and y(x_i) ≈ w_i^.

Then

y(x_i+1) ≈ w_i +

Z xi+1

xi

P (x)dx

s k

= s(s − 1) · · · (s − k + 1) k!

(51)

IVP of ODE

¹⁹

y(xi+1) − y(xi) =

Z x_i+1 xi

y⁰(x)dx =

Z x_i+1 xi

f (x, y(x))dx

⇒ y(x_i+1) = y(xi) +

Z xi+1

xi

f (x, y(x))dx

Z xi+1

xi

P (x)dx

s k

= s(s − 1) · · · (s − k + 1) k!

(52)

IVP of ODE

¹⁹

y(xi+1) − y(xi) =

Z x_i+1 xi

y⁰(x)dx =

Z x_i+1 xi

f (x, y(x))dx

⇒ y(x_i+1) = y(xi) +

Z xi+1

xi

f (x, y(x))dx

Z xi+1

xi

P (x)dx

➳ Adams-Bashforth explicit m-step technique:

The binomial formula is defined as

s k

= s(s − 1) · · · (s − k + 1) k!

(53)

IVP of ODE

¹⁹

y(xi+1) − y(xi) =

Z x_i+1 xi

y⁰(x)dx =

Z x_i+1 xi

f (x, y(x))dx

⇒ y(x_i+1) = y(xi) +

Z xi+1

xi

f (x, y(x))dx

Z xi+1

xi

P (x)dx

➳ Adams-Bashforth explicit m-step technique:

The binomial formula is defined as

s k

= s(s − 1) · · · (s − k + 1) k!

(54)

IVP of ODE

²⁰

We introduce the backward difference notation ∇

∇f (xi) = f (xi) − f (xi−1)

and

∇^kf (x_i) = ∇^k−1f (x_i) − ∇^k−1f (x_i−1) = ∇ ∇^k−1f (x_i) ,

for i = 0, 1, . . . , n − 1^.

Let P_m−1(x) be the Newton backward difference polynomial through (xi, f (xi, y(xi)))^, (x_i−1, f (x_i−1, y(x_i−1)))^, . . .^, (x_i+1−m, f (x_i+1−m, y(x_i+1−m)))^{, i.e.,}

Pm−1(x) =

m−1

X

k=0

(−1)^k





−s k



∇^kf (xi, y(xi))

where x = xi + sh^{. Then}

(55)

IVP of ODE

²⁰

∇f (xi) = f (xi) − f (xi−1)

and

for i = 0, 1, . . . , n − 1^.

Let P_m−1(x) be the Newton backward difference polynomial through (xi, f (xi, y(xi)))^, (x_i−1, f (x_i−1, y(x_i−1)))^, . . .^, (x_i+1−m, f (x_i+1−m, y(x_i+1−m)))^,

i.e.,

Pm−1(x) =

m−1

X

k=0

(−1)^k





−s k



(56)

IVP of ODE

²⁰

∇f (xi) = f (xi) − f (xi−1)

and

for i = 0, 1, . . . , n − 1^.

Let P_m−1(x) be the Newton backward difference polynomial through (xi, f (xi, y(xi)))^, (x_i−1, f (x_i−1, y(x_i−1)))^, . . .^, (x_i+1−m, f (x_i+1−m, y(x_i+1−m)))^{, i.e.,}

Pm−1(x) =

m−1

X

k=0

(−1)^k





−s k

