# Session 8: ODE stability; stiff system¶

Date: 11/06/2017, Monday

```
In [1]:
```

```
format compact
```

## ODE Stability¶

### Problem statement¶

Consider a simple ODE

where \(a\) is a **positive** real number.

It can be solved analytically:

So the analytical solution decays exponentially. A numerical solution is
allowed to have some error, but it should also be decaying. If the
solution is instead growing, we say it is **unstable**.

### Explicit scheme¶

Forward Eulerian scheme for this problem is

The iteration is given by

The solution at the k-th time step can be written out explicitly

We all know that when \(h\) gets smaller the numerical solution will
converge to the true result. But here we are interested in relatively
large \(h\), i.e. we wonder **what’s the largest possible**
\(h\) **we can use while still getting an OK result**.

Here are the solutions with different \(h\).

```
In [2]:
```

```
%plot -s 400,300
% set parameters, can use different value
a = 1;
y0 = 1;
tspan = 10;
% build function
f_true = @(t) y0*exp(-a*t); % true answer
f_forward = @(k,h) (1-h*a).^k * y0; % forward Euler solution
% plot true result
t_ar = linspace(0, tspan); % for plotting
y_true = f_true(t_ar);
plot(t_ar, y_true)
hold on
% plot Forward Euler solution
for h=[2.5, 1.5, 0.5] % try different step size
kmax = round(tspan/h);
k_ar = 0:kmax;
t_ar = k_ar*h;
% we are not doing Forward Euler iteration here
% since the expression is known, we can directly
% get the entire time series
y_forward = f_forward(k_ar, h);
plot(t_ar, y_forward, '-o')
end
% tweak details
ylim([-1,1]);
xlabel('t');ylabel('y');
legend('true solution', 'h=2.5', 'h=1.5', 'h=0.5' , 'Location', 'Best')
```

There all 3 typical regimes, determined by the magnitude of \([1-ha]\) inside the expression \(y_k= (1-ha)^ky_0\).

**Small**\(h\): \(0 \le [1-ha] < 1\)

In this case, \((1-ha)^k\) will decay with \(k\) and always be positve.

We can solve for the range of \(h\):

This corresponds to h=0.5 in the above figure.

**Medium**\(h\): \(-1 \le [1-ha] < 0\)

In this case, the absolute value of \((1-ha)^k\) will decay with \(k\), but it oscillates between negative and postive. This is not desirable, but not that bad, as our solution doesn’t blow up.

We can solve for the range of \(h\):

This corresponds to h=1.5 in the above figure.

**Large**\(h\): \([1-ha] < -1\)

In this case, the absolute value of \((1-ha)^k\) will **increase**
with \(k\). That’s the worst case because the true solution will be
**decaying**, but our numerical solution insteads gives **exponential
growth**. Here our numerical scheme is totally wrong, not just
inaccurate.

We can solve for the range of \(h\):

This corresponds to h=2.5 in the above figure.

The take-away is, **to obtain a stable solution, the time step size**
\(h\) **needs to be small enough. The time step requirement depends
on** \(a\), **i.e. how fast the system is changing.** If \(a\)
is large, i.e. the system is changing rapidly, then \(h\) has to be
small enough (\(h<1/a\) or \(h<2/a\), depends on what your want)
to capture this fast change.

### Implicit scheme¶

But sometimes we really want to use a large step size. For example, we might only care about the steady state where \(t\) is very large, so we would like to quickly jump to the steady state with very few number of iterations. Implicit method allows us to use a large time step while still keep the solution stable.

Backward Eulerian scheme for this problem is

In general, the right hand-side would be some function
\(f(y_{k+1})\), and we need to solve a nonlinear equation to obtain
\(y_{k+1}\). This is why this method is called **implicit**. But in
this problem here we happen to have a linear function, so we can still
write out the iteration explicitly.

The solution at the k-th time step can be written as

Here are the solutions with different \(h\).

```
In [3]:
```

```
%plot -s 400,300
% set parameters, can use different value
a = 1;
y0 = 1;
tspan = 10;
% build function
f_true = @(t) y0*exp(-a*t); % true answer
f_backward = @(k,h) (1+h*a).^(-k) * y0; % backward Euler solution
% plot true result
t_ar = linspace(0, tspan); % for plotting
y_true = f_true(t_ar);
plot(t_ar, y_true)
hold on
% plot Forward Euler solution
for h=[2.5, 1.5, 0.5] % try different step size
kmax = round(tspan/h);
k_ar = 0:kmax;
t_ar = k_ar*h;
% we are not doing Forward Euler iteration here
% since the expression is known, we can directly
% get the entire time series
y_backward = f_backward(k_ar, h);
plot(t_ar, y_backward, '-o')
end
% tweak details
ylim([-0.2,1]);
xlabel('t');ylabel('y');
legend('true solution', 'h=2.5', 'h=1.5', 'h=0.5' , 'Location', 'Best')
```

Since \(\frac{1}{1+ha}\) is always smaller than 1 for any positive
\(h\) and postive \(a\), \(y_k= \frac{y_k}{(1+ha)^k}\) will
always decay. So we don’t have the instability problem as in the
explicit method. A large \(h\) simply gives **inaccurate** results,
but not **terribly wrong** results.

According to no free lunch theorm, implicit methods must have some additional costs (half joking. that’s another theorm). The cost for an implicit method is solving a nonlinear system. In general we will have \(f(y, t)\) on the right-hand side of the ODE, not simply \(-ay\).

### General form¶

Using \(-ay\) on the right-hand side allows a simple analysis, but the idea of ODE stability/instability applies to general ODEs. For example considering a system like

You can use the typical magnitude of \(f(t)\) as the “\(a\)” in the previous analysis.

Even for

You can consider the typical magnitude of \(yf(t)\)

## Stiff system¶

Consider an ODE system

Using an explicit method, the time step requirement for the first equation is \(h<1/a_1\), while the requirement for the second one is \(h<1/a_2\). If \(a_1 >> a_2\), we have to use a quite small \(h\) to accomodate the first requirement, but that’s an overkill for the second equation. You will be using too many unnecessary time steps to solve \(y_2\)

We can define the **stiff ratio** \(r=\frac{a_1}{a_2}\). A system is
very stiff if \(r\) is very large. With explicit methods you often
need an unnecessarily large amount of time steps to ensure stability.
Implicit methods are particularly useful for a stiff system because it
has no instability problem.

You might want to solve two equations separately so we can use a larger time step for the second equation to save computing power. But it’s not that easy because in real examples the two equations are often intertwined: