Numerical Integration of Ordinary Differential Equations
This chapter presents an overview of numerical integration techniques for solving ODE systems, as implemented in Matlab and COMSOL. These techniques are broadly classified into one-step and multistep methods. One-step methods include forward and backward
- PDF / 758,900 Bytes
- 50 Pages / 439.37 x 666.142 pts Page_size
- 60 Downloads / 271 Views
Numerical Integration of Ordinary Differential Equations
Almost all lumped-parameter ODE models encountered in physics, biology or medicine cannot be solved in closed-loop form, but require numerical integration to find their solution. As we have seen from the previous chapter, Matlab provides an extensive suite of ODE solvers including ode15s and ode45, and COMSOL also employs its own ODE solvers based on the backward differential formula (BDF) and generalized-α methods. This chapter provides an overview of numerical methods commonly used for integrating ODE systems, including those algorithms used by Matlab and COMSOL.
3.1 Taylor’s Theorem The basis of all numerical integration algorithms is Taylor’s theorem, which can be stated for a function of a single variable as: Theorem 3.1 Let f be a real-valued function of a single variable with n + 1 continuous derivatives in the closed interval [t, t + h]. Then there exists some ξ ∈ [t, t + h] such that f (t + h) = f (t) + h f (t) +
h 2 h n (n) h n+1 f (t) + · · · + f (t) + f (n+1) (ξ ) 2! n! (n + 1)! (3.1)
© Springer-Verlag Berlin Heidelberg 2017 S. Dokos, Modelling Organs, Tissues, Cells and Devices, Lecture Notes in Bioengineering, DOI 10.1007/978-3-642-54801-7_3
55
56
3 Numerical Integration of Ordinary Differential Equations
Proof We begin with the fundamental theorem of calculus:
b
f (τ ) dτ = f (b) − f (a)
a
Using a = t, b = t + h and re-arranging: R1
t+h
f (t + h) = f (t) +
f (τ ) dτ
t
The integral on the right (remainder term R1 ) can be expanded using the integration by parts formula:
b
b u · v dτ = u v dτ −
a
b
a
a
du dτ
v dτ
dτ
(3.2)
Using u = t + h − τ and v = f (τ ), we have:
t+h t
t+h (t + h − τ ) f (τ ) dτ = (t + h − τ ) f (τ ) t −
t+h
= −h f (t) +
t+h
− f (τ ) dτ
t
f (τ ) dτ
t
Hence, R1 =
t+h
R2
t+h
f (τ ) dτ = h f (t) +
t
(t + h − τ ) f (τ ) dτ
t
We can continue expanding R2 using integration by parts (Eq. 3.2), now with u = (t+h−τ )2 and v = f (τ ): 2!
t+h t
t+h t+h (t + h − τ )2 − −(t + h − τ ) f (τ ) dτ f (τ ) 2! t t t+h h2 = − f (t) + (t + h − τ ) f (τ ) dτ 2! t
(t + h − τ )2 f (τ ) dτ = 2!
Hence, R2 = t
t+h
(t + h − τ ) f (τ ) dτ =
h 2 f (t) + 2!
t
R3
t+h
(t + h − τ )2 f (τ ) dτ 2!
3.1 Taylor’s Theorem
57
We can continue indefinitely expand the remainder terms integration by parts to obtain
f (t + h) = f (t) +
n
i=1
h i (i) f (t) + i!
Rn+1
t+h
t
(t + h − τ )n (n+1) f (τ ) dτ n!
(3.3)
To verify that Eq. 3.3 is true, we make use of the principle of mathematical induction. Assume it is true for some value n = k. That is, f (t + h) = f (t) +
i=1
with
k
hi
t+h
Rk+1 = t
i!
f (i) (t) + Rk+1
(t + h − τ )k (k+1) f (τ ) dτ k!
We then show that Eq. 3.3 also holds for k + 1. This is equivalent to showing that h k+1 f (k+1) (t) + Rk+2 . To do this, we use integration by parts (Eq. 3.2) with Rk+1 = (k+1)! u= t
(t+h−τ )k+1 (k+1)!
Data Loading...
