Biomedical Engineering Reference
In-Depth Information
point halfway across the interval. This is our trial step. Using the derivative
at this estimated point, we advance the dynamics across the full width of the
interval. By this procedure, we obtain second-order accuracy.
We can illustrate this method with a one-dimensional example: the inte-
gration of
x
=
f
(
x
). We define an estimate of the value of the variable at an
instant
t
+∆
t
by
x
estimated
(
t
+∆
t
)=
x
(
t
)+∆
tf
x
+
1
2
∆
tf
(
x
)
,
(7.6)
where
f
(
x
+
2
∆
tf
(
x
)) is
x
evaluated at the point
x
+(1
/
2)∆
tf
(
x
). On the
other hand, we can expand
x
(
t
+∆
t
) as a Taylor series in ∆
t
, and use it to
compute the
error
,
error
=
x
(
t
+∆
t
)
−
x
estimated
(
t
+∆
t
)
.
(7.7)
Expanding
x
(
t
+∆
t
) as a Taylor series in ∆
t
, we obtain
error
=
O
(∆
t
3
),
since the quadratic terms cancel out.
It is natural to ask then why we should not continue with the process of
estimating the time derivatives at yet more points, and obtain in that way an
approximation which is better and better for a given size of time step. There
is nothing wrong with that argument . . . except that, at some point, what
we gain on one side is lost on the other side, since many functions have to
be evaluated. A good compromise is known as the fourth-order Runge-Kutta
method, which, for a system of equations written as
x
=
f
(
x
)
,
(7.8)
estimates
x
n
+1
=
x
n
+
1
6
(
k
1
+2
k
2
+2
k
3
+
k
4
)
,
(7.9)
where
k
1
=
f
(
x
n
)∆
t,
(7.10)
k
2
=
f
x
n
+
1
∆
t,
2
k
1
(7.11)
k
3
=
f
x
n
+
1
∆
t,
2
k
2
(7.12)
k
4
=
f
(
x
n
+
k
3
)∆
t.
(7.13)
Typically, very e
cient routines are used to implement these techniques. The
interested reader is encouraged to consult classical references such as
Numer-
ical Recipes in C
[Press et al. 1992], which contain widely tested subroutines
that compute the quantities above. In this way, writing a code that approx-
imates the solutions of our differential equations is as easy as implementing
a loop in which, at each step, the subroutine is called.
