Ekstra opgaver

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
plt.rc('font',size=16)
from scipy.integrate import solve_ivp

t0 = 0
tmax = 5 * np.pi
y0 = 1

dydt = lambda t,y: 0.9 * y**2 * np.cos(t)

def solve_it(delta_t=0.01):
    t = t0
    y = y0
    ts = [t]
    ys = [y]
    
    for _ in np.arange(t,tmax,delta_t):
       y += dydt(t,y) * delta_t
       t += delta_t
       ts.append(t)
       ys.append(y)

    return ts,ys


fig, ax = plt.subplots()

for delta_t in [0.01, 0.05, 0.25]:
    ts,ys = solve_it(delta_t)
    ax.plot(ts,ys,label='$\Delta t={:6.2f}$'.format(delta_t))

ax.legend()

def solve_it(delta_t=0.01):
    t = t0
    y = y0
    ts = [t]
    ys = [y]
    
    for _ in np.arange(t,tmax,delta_t):
        s1 = dydt(t,y)
        s2 = dydt(t + delta_t/2,y + s1 * delta_t/2)
        
        y += s2 * delta_t
        t += delta_t
        ts.append(t)
        ys.append(y)

    return ts,ys

plt.rc('font',size=16)
fig, ax = plt.subplots()

for delta_t in [0.01, 0.05, 0.25]:
    ts,ys = solve_it(delta_t)
    ax.plot(ts,ys,label='$\Delta t={:6.2f}$'.format(delta_t))

ax.legend()

plt.rc('font',size=16)
fig, ax = plt.subplots()

for ms in [0.01, 0.05, 0.25]:
    mysol = solve_ivp(dydt, [t0,tmax], [y0], max_step=ms)
    ts = mysol.t
    ys = mysol.y[0]
    ax.plot(ts,ys,label='$\mathrm{{ms}}={:6.2f}$'.format(ms))

ax.legend()

def my_ivp_solver(dydt,trange,yinit,max_step=0.01):
    # kald input det samme som før
    t0 = trange[0]
    tmax = trange[1]
    y0 = yinit[0]
    delta_t = max_step

    # gør det samme som før
    t = t0
    y = y0
    ts = [t]
    ys = [y]
    
    for _ in np.arange(t0,tmax,delta_t):
        s1 = dydt(t,y)
        s2 = dydt(t + delta_t/2,y + s1 * delta_t/2)
        
        y += s2 * delta_t
        t += delta_t
        ts.append(t)
        ys.append(y)

    f = lambda x: None
    f.t = np.array(ts)
    f.y = np.array([ys])
    return f

plt.rc('font',size=16)
fig, ax = plt.subplots()

for ms in [0.01, 0.05, 0.25]:
    mysol = my_ivp_solver(dydt, [t0,tmax], [y0], max_step=ms)
    ts = mysol.t
    ys = mysol.y[0]
    ax.plot(ts,ys,label='$\mathrm{{ms}}={:6.2f}$'.format(ms))

ax.legend()