{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# ODE System"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider a system of two coupled ODE\n",
"$$\n",
"y' = Ay, \\qquad A = \\begin{bmatrix}\n",
"-100 & 1 \\\\\n",
"0 & -1/10 \\end{bmatrix}\n",
"$$\n",
"with initial condition\n",
"$$\n",
"y(0) = y_0\n",
"$$"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The $\\theta$ scheme is given by\n",
"$$\n",
"y^{n+1} = y^n + h [(1-\\theta) Ay^n + \\theta Ay^{n+1}]\n",
"$$\n",
"i.e.\n",
"$$\n",
"[I-\\theta hA] y^{n+1} = y^n + (1-\\theta)h A y^n\n",
"$$\n",
"The following function implements the $\\theta$ scheme."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def solve(y0,h,T,theta):\n",
" A = np.matrix([[-100.0, 1.0], [0.0, -0.1]])\n",
" A1 = np.eye(2) - h*theta*A\n",
" A1 = np.linalg.inv(A1)\n",
" A2 = np.eye(2) + (1.0-theta)*h*A\n",
" M = A1.dot(A2)\n",
" N = int(T/h)\n",
" y = np.ndarray((N,2))\n",
" t = np.zeros(N)\n",
" y[0,:]=y0\n",
" t[0] = 0.0\n",
" for i in range(N-1):\n",
" y[i+1,:] = M.dot(y[i,:])\n",
" t[i+1] = t[i] + h\n",
" plt.figure(figsize=(9,4))\n",
" plt.subplot(1,2,1)\n",
" plt.plot(t,y[:,0])\n",
" plt.xlabel('t')\n",
" plt.title('First component')\n",
" plt.subplot(1,2,2)\n",
" plt.plot(t,y[:,1])\n",
" plt.xlabel('t')\n",
" plt.title('Second component')\n",
" \n",
"# Initial condition\n",
"y0 = np.array([10.0/999.0,1.0])\n",
"\n",
"# Final time\n",
"T = 25.0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward Euler, $\\theta=0$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The stability condition is\n",
"$$\n",
"|1-100h| < 1, \\qquad |1-h/10| < 1\n",
"$$\n",
"which means that we need to choose $h < 1/50$. We first try with $h=0.01$"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.01\n",
"solve(y0,h,T,0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We now try with $h=0.021$ which is slightly larger than the stability limit."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.021\n",
"solve(y0,h,T,0.0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Trapezoidal scheme, $\\theta = 0.5$"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.021\n",
"solve(y0,h,T,0.5)"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 2
}