- Form of ODE: $$\frac{dy}{dt} = f(y,t)$$
- (First order ODE)
- Order of an ODE is the highest derivative present.
- Examples of first, second, third order ODEs?
- Order of the ODE is different than the $\mathcal{O}$ order of a truncation error.
- Boundary conditions (BCs):
- Need the same number of boundary conditions as the order of the ODE.
- For $y=y(t)$, we call the BC an initial condition (I.C.).
- Most first order ODEs use initial conditions.
- Doesn’t have to be time
- The solution $y(t)$ satisfies both the ODE and the BCs.
- Linear ODE: $y$ and all derivatives of $y$ are to the $1^{st}$ power, no products of y with itself or derivatives.
- $y^{\prime} + \alpha y = f(t)$
- linear for $\alpha(t)$, but not for $\alpha(y,t)$.
- $yy^{\prime} + \alpha y = 0$
- nonlinear
- $y^{\prime} + \alpha t y = f(t)$
- linear
- $y^{\prime} + \alpha y = f(t)$
- System of ODEs: Multiple dependent variables, one independent variable.
- $dw/dt = f(w,z,t)$,
- $dz/dt = g(w,z,t)$.
- We would normally write these in vector notation: $$\frac{d\vec{y}}{dt} = \vec{f}(\vec{y},t),$$
- Or, just $$\frac{dy}{dt} = f(y,t),$$ where it’s understood that we are dealing with vector quantities.
- Here, $f$ is a function that takes a vector of variables and returns a vector of rates of those variables.
- 2 classes
- Initial value problems: marching methods, usually time, usually first order ODE.
- Boundary value problems:marching or equilibrium methods, usually higher order ODE (2 end points, solutions on spatial domains).
- Physical classes
- Propagation: IVP, order $\ge 1$.
- Equilibrium: BVP, order $\ge 2$.
- Eigenproblems: solutions only for special parameters.
Methods
- Finite difference methods are emphasized here
- Others include the class of Galerkin methods (finite element…)
- Transform a continuous calculus problem into an algebraic problem.
- Discretize the independent variable into a grid
- Approximate derivatives with F.D. approximations (Taylor series).
- Solve
IVP: focus on single point methods (advance using only the previous point), as opposed to multi-step methods.
Explicit Euler method
$$\frac{dy}{dt} = f(y,t),$$ $$y(t=0) = y_0.$$
- Use a first order forward difference approximation to the derivative $$\frac{y_{k+1} - y_k}{\Delta t} = f(y_k,t),$$
-
First order forward difference.
-
Explicit since we can solve for the new point explicitly in terms of the old point.
-
Taylor Series: $$ y_{k+1} = y_k + y_k^{\prime}\Delta t + \frac{1}{2}y_k^{\prime\prime}\Delta t^2 + \ldots$$
-
The error per step is $\mathcal{O}(\Delta t^2).$
-
The cumulative error is $\mathcal{O}(\Delta t).$
- If we take $n=t_{end}/\Delta t$ steps, the total error is $k\mathcal{O}(\Delta t^2)=\frac{t_{end}}{\Delta t}\mathcal{O}(\Delta t^2) = \mathcal{O}(\Delta t).$
-
Generally, the global error is the order of the F.D. approximation to the derivative.
-
So we say that the explicit Euler (EE) method is globally first order accurate.
- This means that the error is proportional to our step size.
-
EE method is very easy to implement in Excel.
Example 1 (EE)
$$\frac{dy}{dt} = -ay,$$ $$y_0 = 1,$$ $$a=1.1.$$ This equation has the exact solution $$y(t) = y_0e^{-at}.$$
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.optimize import fsolve
def ode_EE(f, y0, t):
ns = len(t)-1
y = np.empty(ns+1)
y[0] = y0
for k in range(ns):
Δt = t[k+1] - t[k]
y[k+1] = y[k] + Δt*f(y[k],t[k])
return y
def myfunc(y,t):
return -1.1*y
#---------------
y0 = 1
tend = 6.0
t = np.linspace(0, tend, 20)
y = ode_EE(myfunc, y0, t)
plt.rc("font", size=14)
plt.plot(t,y)
plt.plot(t, y0*np.exp(-1.1*t), 'r.')
plt.legend(['EE', 'Exact'], frameon=False)
plt.xlabel('t')
plt.ylabel('y');
Implicit Euler method
- Instead of approximating point $y_{k+1}$ in terms of point $y_k$, approximate point $y_k$ in terms of point $y_{k+1}$. $$y_k = y_{k+1} - y_{k+1}^{\prime} \Delta t + \frac{1}{2}y_{k+1}^{\prime\prime}\Delta t^2 + \ldots$$
- This can be rearranged to the following, using $y^{\prime}{k+1} = f(y{k+1},t_{k+1}).$
- This is like EE, except instead of evaluating the rate at point $k$, it is evaluated at point $k+1$.
- This makes the method implicit since the equation for $y_{k+1}$ is implicit in $y_{k+1}$.
- For linear ODEs, the Implicit Euler (IE) method can be solved directly for $y_{k+1}$.
- For nonlinear ODEs, each advancement step requires solution of a nonlinear system of the form $F(y_{k+1}) = 0 = y_{k+1}-y_{k}-\Delta tf(y_{k+1},t_{k+1}).$
- Globally first order accurate: $\mathcal{O}(\Delta t)$.
Example 2 (IE)
- Solve the same problem as we did before: $dy/dt = -ay$ with $a=1.1$ and $y_0=1$.
- Solve the blue IE iteration equation for $y_{k+1}$: $$ y_{k+1} = \frac{y_k}{1+a\Delta t}.$$
y0 = 1
a = 1.1
tend = 6.0
t = np.linspace(0.0, tend, 20)
ns = len(t) - 1
y = np.zeros(ns+1)
y[0] = y0
for k in range(ns): # have n points, take n-1 steps
Δt = t[k+1]-t[k]
y[k+1] = y[k]/(1+a*Δt)
plt.rc("font", size=14)
plt.plot(t,y)
plt.plot(t, y0*np.exp(-a*t), 'r.')
plt.legend(['IE', 'Exact'], frameon=False)
plt.xlabel('t')
plt.ylabel('y')
Text(0, 0.5, 'y')
Example 3 (EE, IE nonlinear)
- Solve the following nonlinear ODE using EE and IE $$\frac{dy}{dt} = -y^{0.8},$$ $$y_0 = 1.$$
def ode_IE(f, y0, t):
ns = len(t) - 1
y = np.empty(ns+1)
y[0] = y0
def F(ykp1):
return ykp1-y[k]-Δt*f(ykp1,t[k+1])
for k in range(ns):
Δt = t[k+1]-t[k]
y[k+1] = fsolve(F, y[k])[0]
return y
def myfunc(y,t):
return -y**0.8
#------------------
y0 = 1.0
tend = 3.0
t = np.linspace(0.0, tend, 20)
n = len(t)
yEE = ode_EE(myfunc,y0,t)
yIE = ode_IE(myfunc,y0,t)
#---------------------
plt.rc("font", size=14)
plt.plot(t,yEE,'b-')
plt.plot(t,yIE,'g-')
plt.plot(t, (1-0.2*t)**(1/0.2), 'r.')
plt.legend(['EE', 'IE', 'Exact'], frameon=False)
plt.xlabel('t')
plt.ylabel('y');
