# ODE

We will solve the following linear Cauchy model

\begin{align}
y^{\prime}(t) &= \lambda y(t)\\
y(0) & = 1
\end{align}

whose exact solution is

$$
y(t) = e^{\lambda t}
$$



```python
%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import scipy.linalg
import numpy.linalg
```


```python
l = -5.   # lamba (see formula above)
t0 = 0.   # initial time
tf = 10.  # final time
y0 = 1.   # boundary condition

s = linspace(t0,tf,5000)
exact = lambda x: exp(l*x)  # exact solution
```

### Forward Euler

$$
\frac{y_{n}-y_{n-1}}{h} = f(y_{n-1}, t_{n-1})
$$


```python
def fe(l,y0,t0,tf,h):
    timesteps = arange(t0,tf+1e-10, h)
    sol = zeros_like(timesteps)
    sol[0] = y0
    for i in range(1,len(sol)):
        sol[i] = sol[i-1]*(1+l*h)
    
    return sol, timesteps

y, t = fe(l,y0,t0,tf,0.1)
_ = plot(t,y, 'o-')
_ = plot(s,exact(s))

error = numpy.linalg.norm(exact(t) - y, 2)
print(error)
```

### Backward Euler

$$
\frac{y_{n}-y_{n-1}}{h} = f(y_{n}, t_{n})
$$


```python
def be(l,y0,t0,tf,h):
    timesteps = arange(t0, tf+1e-10, h)
    sol = zeros_like(timesteps)
    sol[0] = y0
    for i in range(1, len(sol)):
        sol[i] = sol[i-1]/(1-l*h)

    return sol, timesteps

y,t = be(l, y0, t0, tf, 0.1)
_ = plot(t, y, 'o-')
_ = plot(s, exact(s))

error = numpy.linalg.norm(exact(t)-y, infty)
print(error)
```

### $\theta$-method

$$
\frac{y_{n}-y_{n-1}}{h} = \theta\, f(y_{n}, t_{n}) + (1-\theta)\,f(y_{n-1}, t_{n-1})
$$


```python
def tm(theta,l,y0,t0,tf,h):
    timesteps = arange(t0, tf+1e-10, h)
    sol = zeros_like(timesteps)
    sol[0] = y0
    for i in range(1, len(sol)):
        sol[i] = theta*sol[i-1]/(1-l*h) + (1. - theta)*sol[i-1]*(1+l*h)

    return sol, timesteps

y,t = tm(0.5, l, y0, t0, tf, 0.1)
_ = plot(t, y, 'o-')
_ = plot(s, exact(s))

error = numpy.linalg.norm(exact(t) - y, infty)
print(error)
```

### Simple adaptive time stepper

For each time step:
- Compute solution with CN
- Compute solution with BE
- Check the difference
- If the difference satisfy a given tolerance:
    - keep the solution of higher order
    - double the step size
    - go to the next step
- Else:
    - half the step size and repeat the time step


```python
def adaptive(l,y0,t0,tf,h0, hmax=0.9,tol=1e-3):
    sol = []
    sol.append(y0)
    t = []
    t.append(t0)
    h = h0
    while t[-1] < tf:
        #print 'current t =', t[-1], '            h=', h
        current_sol = sol[-1]
        current_t = t[-1]
        sol_cn, _ = tm(0.5,l,current_sol,current_t, current_t + h, h)
        sol_be, _ = tm(1.,l,current_sol,current_t, current_t + h, h)
        
        if (abs(sol_cn[-1] - sol_be[-1]) < tol): #accept
            sol.append(sol_cn[-1])
            t.append(current_t+h)
            h *= 2.
            if h > hmax:
                h=hmax
        else:
            h /= 2.
            
    return sol, t

y,t = adaptive(l,y0,t0,tf,0.9)
_ = plot(t,y, 'o-')
_ = plot(s,exact(array(s)))

error = numpy.linalg.norm(exact(array(t)) - y, infty)
print(error, len(y))
```

40ee35fda98d92afd1479d68c4221de22a8c4e31
ipynb
Jupyter Notebook
notebooks/GR_lab06a_ODEs_1.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
notebooks/GR_lab06a_ODEs_1.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
notebooks/GR_lab06a_ODEs_1.ipynb
mapenzo-ph/numerical-analysis-2021-2022
952808d7fa7a9e718274592104be24882acef3ec
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES