Numerical solution methods - Essential Mathematics for Chemical Engineers

Some ODEs are hard (or impossible!) to solve analytically... or maybe you just don’t want to! Numerical approximations to the solution can be very valuable in these cases.

For example, let’s solve the initial value problem:

y' - y = t, \qquad y(0) = 0

(3)

using Euler’s method with a step size $\Delta t = 0.2$ . First, rewrite in explicit form:

y' = t + y = f(t, y)

(4)

Then, form a table and update

$n$	$t$	$y$	$f$
0	0	0	0
1	0.2	0	0.2
2	0.4	0.04	0.44
3	0.6	0.128	0.728
4	0.8	0.274	1.074
5	1.0	0.489

Let’s compare our numerical result to the analytical solution to the same ODE. This solution can be obtained using an integrating factor:

F = e^{\int -1 \d{t}} = e^{-t}

(5)

Integrate by parts:

\begin{align} \int F r \d{t} &= \int t e^{-t} \d{t} \\ &= -t e^{-t} + \int e^{-t} \d{t} \\ &= -(t + 1) e^{-t} \end{align}

(6)

So the general solution is:

y(t) = e^t \left[-(t + 1) e^{-t} + c \right] = c e^x - t - 1

(7)

Solve for c using the initial condition:

y(0) = c - 1 = 0 \to c = 1

(8)

y(t) = e^t - t - 1

(9)

Comparing the numerical and analytical solutions at the same time points by computing the absolute error between them gives:

$t$	Numerical	Analytical	Error
0	0	0	0
0.2	0	0.021	0.021
0.4	0.04	0.092	0.052
0.6	0.128	0.222	0.094
0.8	0.274	0.426	0.152
1.0	0.489	0.718	0.229

The error increases as more timesteps are run due to accumulation of error. A smaller timestep typically produces a better solution, but at the cost of doing additional calculations. We will discuss next how the error depends on the step size.

Error in approximation¶

Because we are using truncated Taylor series, we accrue an error:

$Euler's method error$

The error in a single step is $O(\Delta t^2)$ . This is called the local error. To cover a finite time t, $t/\Delta t$ steps are required so the global (total) error is $O(\Delta t)$ . To improve accuracy, you should decrease $\Delta t$ , but this takes more work. Some $\Delta t$ will also totally fail! This is the field of numerical stability (not covered here).

Skill builder problems¶

Solution to Exercise 1

The ODE is already in explicit form, so $f(t, y) = y + 5 \sin(2\pi t)$ .

Euler’s method with $\Delta t = 0.1$ :

$n$	$t$	$y$	$f$
0	0.0	1.000	1.000
1	0.1	1.100	4.039
2	0.2	1.504	6.257
3	0.3	2.130	6.895
4	0.4	2.818	5.757
5	0.5	3.394	3.394
6	0.6	3.733	0.795
7	0.7	3.813	-0.942
8	0.8	3.719	-1.037
9	0.9	3.615	0.676
10	1.0	3.683

gives $y(1) \approx 3.683$ .

Euler’s method with $\Delta t = 0.2$ :

$n$	$t$	$y$	$f$
0	0.0	1.000	1.000
1	0.2	1.200	6.955
2	0.4	2.591	5.330
3	0.6	3.657	0.576
4	0.8	3.772	-1.155
5	1.0	3.322

gives $y(1) = 3.322$ .

To find the analytical solution, note that this is a linear first-order ODE

y' - y = 5 \sin(2\pi t)

(10)

with $p = -1$ and $r = 5 \sin (2\pi t)$ . The integrating factor and required integral are:

\begin{align} F &= e^{\int p \d{t}} = e^{-t} \\ \int F r \d{t} &= 5 \int e^{-t} \sin(2\pi t) \d{t} \\ &= \frac{5 e^{-t}}{4\pi^2 + 1}\left[ -\sin(2\pi t) - 2\pi \cos(2\pi t)\right] \end{align}

(11)

where the integral was evaluated using a table of integrals. Then, solve for y

\begin{align} y &= \frac{1}{F} \left[\int F r \d{t} + c \right] \\ &= -\frac{5}{4\pi^2 + 1} \left(\sin(2\pi t) + 2\pi \cos(2\pi t)\right) + c e^{t} \end{align}

(12)

Finally, apply the initial condition:

y(0) = \frac{-10 \pi}{4\pi^2 + 1} + c = 1

(13)

y = \left(1 + \frac{10\pi}{4\pi^2 + 1}\right) e^{t} - \frac{5}{4\pi^2 + 1} \left[\sin(2\pi t) + 2\pi \cos(2\pi t)\right]

(14)

Evaluating this at $t = 1$ gives $y(1) = 4.052$ .

The absolute error is 0.369 for $\Delta t = 0.1$ and 0.730 for $\Delta t = 0.2$ . Hence the error decreases by about a factor of 2 when $\Delta t$ is also decreased by a factor of 2.

5.6 Numerical solution methods

Error in approximation¶

Skill builder problems¶