Initial Value Non-Homogeneous Linear ODE Example
Contents
Initial Value Non-Homogeneous Linear ODE Example¶
Randall Romero Aguilar, PhD
This demo is based on the original Matlab demo accompanying the Computational Economics and Finance 2001 textbook by Mario Miranda and Paul Fackler.
Original (Matlab) CompEcon file: demode03.m
Running this file requires the Python version of CompEcon. This can be installed with pip by running
!pip install compecon --upgrade
Last updated: 2021-Oct-01
About¶
Solve
\[\begin{align*}
\dot{x_1} &= x_1 +12x_2 - 60\\
\dot{x_2} &= -x_1 - 6x_2 + 36
\end{align*}\]
subject to
\[\begin{align*}
x_1(0) &=5\\
x_2(1) &=2\\
t &\in [0, 3]
\end{align*}\]
FORMULATION¶
from compecon import jacobian, ODE
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
Phase Diagram Title and Axis Labels¶
xlabels = ['$x_1$', '$x_2$']
Velocity Function¶
A = np.array([[1, 12],[-1, -6]])
b = np.array([-60,36])
def f(x, A, b):
return ((A @ x).T + b).T
Closed-Form Solution¶
def X(t):
return pd.DataFrame({xlabels[0]: 12 - 48*np.exp(-2*t) + 42*np.exp(-3*t),
xlabels[1]: 4 + 12*np.exp(-2*t) - 14*np.exp(-3*t)},
index = t)
Initial State¶
xinit = np.array([6,2])
Time Horizon¶
T = 3
SOLVE ODE ANALYTICALLY¶
# Time Discretization
N = 100 # number of time nodes
t = np.linspace(0,T, N) # time nodes
# Plot Closed-Form Solution in Time Domain
fig, ax= plt.subplots()
X(t).plot(ax=ax)
ax.set(title='Solution', xlabel='Time')
[Text(0.5, 1.0, 'Solution'), Text(0.5, 0, 'Time')]
SOLVE ODE USING RUNGE-KUTTA METHOD (RK4)¶
# Solve ODE
problem = ODE(f, T, xinit, A, b, labels=xlabels)
problem.rk4(n=N)
problem.x.plot()
<AxesSubplot:>
# Plot Runge-Kutta Approximation Errors
fig, ax = plt.subplots()
(problem.x - X(problem.x.index)).plot(ax=ax)
ax.axhline(0,color='black', ls=':')
ax.set(title='Runge Kutta Approximation Errors', xlabel='Time', ylabel='Error')
[Text(0.5, 1.0, 'Runge Kutta Approximation Errors'),
Text(0.5, 0, 'Time'),
Text(0, 0.5, 'Error')]
problem.x.plot()
<AxesSubplot:>
SOLVE ODE USING COLLOCATION METHOD (ODECOL)¶
# Solve ODE
n = 20; # number of basis functions
problem.solve_collocation(n=n)
# Plot Collocation Approximation Errors
fig, ax = plt.subplots()
t = problem.x.index
(problem.x - X(t)).plot(ax=ax)
ax.axhline(0,color='black', ls=':')
ax.set(title='Collocation Approximation Errors',
xlabel='Time',
ylabel='Error')
[Text(0.5, 1.0, 'Collocation Approximation Errors'),
Text(0.5, 0, 'Time'),
Text(0, 0.5, 'Error')]
STEADY-STATE¶
# Compute Steady State
xstst = -np.linalg.solve(A,b)
print('Steady State')
print(xstst)
print('Eigenvalues')
print(np.linalg.eigvals(jacobian(f,xstst, A,b)))
Steady State
[12. 4.]
Eigenvalues
[-2. -3.]
PHASE DIAGRAM¶
# Plotting Limits
x1lim = [0, 15] # x1 plotting limits
x2lim = [0, 8] # x2 plotting limits
# Compute Separatrix
problem.spx(x0=xstst,T=12)
# Compute Nullclines
x1 = np.linspace(*x1lim, 100)
xnulls = pd.DataFrame({'$\dot{x}_1 =0$': -(A[0,0]*x1 + b[0])/A[0,1],
'$\dot{x}_2 =0$': -(A[1,0]*x1 + b[1])/A[1,1]},
index = x1)
problem.phase(x1lim, x2lim,
xnulls=xnulls,
title='ODE Phase Diagram',
animated=2.5 # requires ffmpeg for video, change to 0 for static image
)
problem.spx(x0=xstst,T=12)