Precept 12: Initial value problems

Today we will consider explicit and implicit methods. In particular, we will consider two explicit methods: Euler's method $$ y_{j+1} = y_j + h f(x_j,y_j) $$ the improved Euler method: $$ y_{j+1} = y_j + \frac{h}{2} \left( f(x_j,y_j) + f(x_{j+1},y_j + h f(x_j, y_j) \right) $$

and two implicit methods: implicit Euler's $$ y_{j+1} = y_j + h f(x_{j+1},y_{j+1}), $$ and the trapezoidal rule method $$ y_{j+1} = y_j + \frac{h}{2} \left( f(x_{j},y_{j}) + f(x_{j+1},y_{j+1}) \right) $$ on two sets of equations.

Task 1:

First, a large non-stiff system of equations. Suppose that $y \in \mathbb{R}^n$ where $n = 500$. Consider the system of differential equations $$ y' = A y $$ where $A \in \mathbb{R}^{n \times n}$ is a symmetric matrix whose eigenvalues are contained in the interval $[1/2,1]$ with the initial condition $$ y(0) = b, $$ where $b \in \mathbb{R}^n$. The matrix $A$ and vector $b$ should be loaded from the .mat file precept12.mat using load (In Python you can use scipy.io.loadmat). Determine $y(1)$ using Euler's method to relative error less than $10^{-4}$. The exact solution is given below.

import scipy.io
import numpy as np
data = scipy.io.loadmat('precept12.mat')
A = data['A']
b = data['b'][:,0]
f = lambda x,y : A @ y

# Exact solution is given by eigenvalue decomposition
[S,U] = np.linalg.eig(A)
y_1_exact = U @ ( np.exp(S) * np.linalg.solve(U,b)) # Exact solution at x =1

# Implement Euler's metohd
def euler():
    return

y = euler(f, 10000, 0, b)
# Compute relative error. Make sure it is less than 1e-4
print('error:')

Task 2

A small (1-dimensional) stiff system of equations. Suppose that $$ y' = 100(\sin x - y), \qquad y(0) = 0. $$ The exact solution is $$ y(x) = \frac{\sin x - 0.01 \cos x + 0.01 e^{-100 x}}{1.0001}. $$ Determine $y(1)$ using implicit Euler with relative accuracy $1e-5$. (You can write down an explicit formula in this situation, no need for Newton's method)

f2 = lambda x,y : 100* ( np.sin(x) - y)
y_exact = lambda x : (np.sin(x) - 0.01 * np.cos(x) + 0.01 * np.exp(-100*x))/1.0001
# Solve with implicit Euler's method to $1e-5$ accuracy.

print('error:', )

Task 3:

Stiff differential equations often have the property that solution curves will converge to a common curve regardless of the initial condition. To visualize this, solve the stiff differential equation from Task 2 for $100$ different values of $y(0)$ between $-1$ and $1$ use $\text{plot}(x,y); \text{hold all}$ in a loop to plot all of these all on the same plot.

import matplotlib.pyplot as plt

# Make plots

Bonus tasks

Use the second the Trapezoidal rule method and Improved Euler method to solve tasks 1 $10^{-8}$ and 2 to $10^{-12}$