We now turn our attention to the numerical solution of ordinary differential equations. We will focus on *linear* ODEs to begin with. The reason we do this is it allows us to view the numerical solution as reducing infinite-dimensional linear equations (the ODE) to finite-dimensional linear systems (matrices).

To emphasize the difference between linear and nonlinear, consider the following two examples:

- Example of a linear ODE:
$$y'(t) - a(t) y(t) =0, y(0)=1$$
Write this as $Ly = 0, y(0)=1$ where $L$ is a
*linear operator*: $$L = D - a(t)$$ where $D$ is the*derivative operator*. This is called a linear operator because, for constant $c_1$ and $c_2$, we have $$L(c_1 y_1 + c_2 y_2) = c_1 L y_1 + c_2 L y_2 = c_1 (y_1' -a(t) y_1) + c_2(y_2' -a (t) y_2)$$

- Example of a nonlinear ODE: the following logistic equation is a nonlinear ODE, since it depends quadratically on $y$:

When the equation is an *initial value problem*, we can solve both linear and nonlinear ODEs using `ode45`

. Here is an example of a linear ODE:

In [1]:

```
using ODE, PyPlot
a=t->-t^2
t,y=ode45((t,y)->a(t)*y,1.,0.:0.1:10.)
plot(t,y);
```

Similarly, we can solve a nonlinear ODE:

In [3]:

```
t,y=ode45((t,y)->y*(1-y),0.5,0.:0.1:10.)
plot(t,y);
```

The distinction between linear and nonlinear becomes greater when we consider other types of problems. The following is a linear intial value problem, since the constraints are specified at a single point, $t=0$:

$$u'' -a(t)u = 0, u(0) = 1, u'(0)=0$$While we can solve this with `ode45`

by rewriting as a system, we can also solve it in ApproxFun:

In [14]:

```
using ApproxFun
D=Derivative() # the derivative operator
B=ivp() # conditions corresponding to u(0),u(1)
a=Fun(t->-t^2,[0.,10.])
u=[B;D^2-a]\[1.,0.]
ApproxFun.plot(u)
```

Out[14]:

However, many problems are not initial value problems, but rather *boundary value problem*. One can imagine this coming up if, say, you observe a bullet at two different points, and know it solves an ODE, and want to recover the solution.

As a simple example, we consider

$$u'' -a(t)u = 0, u(0) = 1, u(10)=0.5$$Because it's linear, we can solve BVPs just as easily as initial value problems:

In [17]:

```
D=Derivative()
B=dirichlet()
a=Fun(t->t^2,[0.,10.])
u=[B;D^2-a]\[1.,0.5]
ApproxFun.plot(u)
```

Out[17]:

Our starting point is to construct *Euler's method* for the ODE

The motivation is the observation that

$$y(h) = y(0) + hy'(0) + {h^2 \over 2} y''(\chi) = y(0) + h y'(0) + O(h^2) = y(0) + h f(0,y(0)) + O(h^2)$$where $h$ is a small step. Thus we have $y(h) \approx y(0) + h f(0,y(0))$, and more generally, $y(t_{k+1}) \approx y(t_k) + h f(t_k,y(t_k))$ for $t_k \triangleq k*h$. We use this to construct an approximation $w_k \approx y(t_k)$ via:

$$w_{k+1} = w_k + h f(t_k,w_k)$$For the linear ODE $$y'(t) = a(t) y(t)$$, this becomes

$$w_0 = y(0), w_{k+1} = (1 + h a(t_k)) w_k$$We observe that this does indeed approximate the solution, by comparing to `ode45`

:

In [21]:

```
a=t->-t^2
# step 1: choose a time step
h=0.1
n=round(Int,8/h)
w=zeros(n+1) # a vector that approximates y at each time step
t=0.:h:10.
w[1]=1.
M=length(w)-1
for k=1:M
w[k+1]=(1 + h*a(t[k]))*w[k]
end
plot(t[1:M+1],w[1:M+1])
t,y=ode45((t,y)->a(t)*y,1.,t)
plot(t,y);
```

Taking a small step size and we get an accurate approximation, though it eventually blows up.

The key observation we make now is that we are in fact solving a lower triangular linear system with forward elimination:

$$ \begin{pmatrix} 1 \cr 1+ha(t_0) & -1 \cr & 1+ha(t_1) & -1 \cr && 1+ha(t_2) & -1 \cr &&&\ddots & \ddots \end{pmatrix} \begin{pmatrix} w_0 \cr w_1 \cr w_2 \cr w_3 \cr \vdots \end{pmatrix} = \begin{pmatrix} 1 \cr 0 \cr 0 \cr 0 \cr \vdots \end{pmatrix} $$