<jupyter_start><jupyter_text>Non-IVP Non-Homogeneous Linear ODE Example**<NAME>, PhD**This demo is based on the original Matlab demo accompanying the  Computational Economics and Finance 2001 textbook by <NAME> and <NAME>.Original (Matlab) CompEcon file: **demode04.m**Running this file requires the Python version of CompEcon. This can be installed with pip by running    !pip install compecon --upgradeLast updated: 2021-Oct-01  AboutSolve\begin{align*}\dot{x_1} &= -1*x_1 - 0.5x_2 + 2\\\dot{x_2} &= -0.5x_2 + 1\end{align*} subject to\begin{align*}x_1(0) &=1\\x_2(1) &=1\\  t &\in [0, 10]\end{align*}  FORMULATION<jupyter_code>from compecon import jacobian, ODE

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt<jupyter_output><empty_output><jupyter_text>Axis Labels<jupyter_code>xlabels  = ['$x_1$','$x_2$']<jupyter_output><empty_output><jupyter_text>Velocity Function<jupyter_code>A = np.array([[-1, -0.5], [0, -0.5]])
b = np.array([2, 1])

def f(x, A, b):
    return ((A @ x).T + b).T<jupyter_output><empty_output><jupyter_text>Boundary Conditions<jupyter_code>bx = np.array([0,1])     # boundary variables
bt = np.array([0,1])     # boundary times
bv = np.array([1,1])     # boundary values<jupyter_output><empty_output><jupyter_text>Time Horizon<jupyter_code>T = 10<jupyter_output><empty_output><jupyter_text>Closed-Form Solution<jupyter_code>def X(t):
    values = np.array([1 - np.exp(1/2-t) + np.exp((1-t)/2), 2-np.exp((1-t)/2)])
    return pd.DataFrame(values.T, index=t, columns=xlabels)<jupyter_output><empty_output><jupyter_text>SOLVE ODE ANALYTICALLY<jupyter_code># Time Discretization
N = 200                 # 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')<jupyter_output><empty_output><jupyter_text>SOLVE ODE USING COLLOCATION METHOD<jupyter_code># Solve ODE
n = 15    # number of basis functions

problem = ODE(f, T, bv, A, b, labels=xlabels)
problem.solve_collocation(n=n, bt=bt, bx=bx, nf=10)

problem.x.plot()
# Plot Collocation Approximation Errors
fig, ax = plt.subplots()

(problem.x - X(problem.x.index)).plot(ax=ax)
ax.axhline(0, color='black', ls=':')
ax.set(title='Collocation Approximation Errors',
       xlabel='Time',
       ylabel='Error');
# Plot Residuals
fig, axs= plt.subplots(1,2)
problem.resid.plot(ax=axs, subplots=True, legend=False)

for ax in axs:
    ax.axhline(0, color='black', ls=':')
    ax.set(xlabel='Time', ylabel='Residual')

    
fig.suptitle('Collocation Residuals')<jupyter_output><empty_output><jupyter_text>STEADY-STATE<jupyter_code># Compute Steady State
xstst = - np.linalg.solve(A,b)
print('Steady State')
print(xstst)
print('Eigenvalues')
print(np.linalg.eigvals(jacobian(f, np.atleast_2d(xstst).T, A, b)))<jupyter_output>Steady State
[1. 2.]
Eigenvalues
[-1.  -0.5]<jupyter_text>PHASE DIAGRAM<jupyter_code># Plotting Limits
x1lim = [0, 2]  # x1 plotting limits
x2lim = [0, 4]  # x2 plotting limits

# Compute Nullclines
x1 = np.linspace(*x1lim, 100)
xnulls = pd.DataFrame({'$x_0$ Nullcline': -(A[0,0]*x1 + b[0])/A[0,1],
                       '$x_1$ Nullcline': -(A[1,0]*x1 + b[1])/A[1,1]},
                     index = x1)

# Plot Phase Diagram
problem.phase(x1lim, x2lim, 
         xnulls=xnulls, 
         xstst=xstst,
         title='ODE Phase Diagram'
        )<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16091281451141018, 'fraction_numerical': 0.031304856641310705, 'mean_word_length': 4.794915254237288, 'pattern_counts': {'":': 0, '<': 43, '<?xml version=': 0, '>': 43, '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': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '205499', 'license': 'MIT', 'n_tokens_mistral': 1456, 'n_tokens_neox': 1379, 'n_words': 273, 'path': '_build/html/_sources/notebooks/ode/04 Non-IVP Non-Homogeneous Linear ODE Example.ipynb', 'repo_name': 'randall-romero/CompEcon-python', 'size': '3447'}