<jupyter_start><jupyter_text>ODEWe 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}$$<jupyter_code>%matplotlib inline
from numpy import *
from matplotlib.pyplot import *
import scipy.linalg
import numpy.linalg

l = -5.
t0 = 0.
tf = 10.
y0 = 1.

s = linspace(t0,tf,5000)

exact = lambda x: exp(l*x)<jupyter_output><empty_output><jupyter_text>Forward Euler$$\frac{y_{n}-y_{n-1}}{h} = f(y_{n-1}, t_{n-1})$$<jupyter_code>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)<jupyter_output>0.211605395525<jupyter_text>Backward Euler$$\frac{y_{n}-y_{n-1}}{h} = f(y_{n}, t_{n})$$<jupyter_code>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)<jupyter_output>0.076565003273<jupyter_text>$\theta$-method$$\frac{y_{n}-y_{n-1}}{h} = \theta\, f(y_{n}, t_{n}) + (1-\theta)\,f(y_{n-1}, t_{n-1})$$<jupyter_code>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)<jupyter_output>0.0276016633937<jupyter_text>Simple adaptive time stepperFor 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<jupyter_code>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)<jupyter_output>0.000817298421905 74
{'fraction_non_alphanumeric': 0.1904916965157929, 'fraction_numerical': 0.051123412569195704, 'mean_word_length': 3.096, 'pattern_counts': {'":': 0, '<': 19, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '5', 'dir': 'jupyter-structured-clean-dedup', 'id': '146791', 'license': 'CC-BY-4.0', 'n_tokens_mistral': 1456, 'n_tokens_neox': 1376, 'n_words': 249, 'path': 'notes/07a_ode.ipynb', 'repo_name': 'lbonaldo/P1.4_seed', 'size': '3071'}