{
"metadata": {
"name": "",
"signature": "sha256:3253b61671a773af8b4435e827d8cc9fee0aefbb056bc3af1ab8f73e37a7bbb6"
},
"nbformat": 3,
"nbformat_minor": 0,
"worksheets": [
{
"cells": [
{
"cell_type": "heading",
"level": 1,
"metadata": {},
"source": [
"ODE with periodic solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Consider the ODE system\n",
"$$\n",
"x' = -y, \\qquad y' = x\n",
"$$\n",
"with initial condition\n",
"$$\n",
"x(0) = 1, \\qquad y(0) = 0\n",
"$$\n",
"The exact solution is\n",
"$$\n",
"x(t) = \\cos(t), \\qquad y(t) = \\sin(t)\n",
"$$\n",
"This solution is periodic. It also has a quadratic invariant\n",
"$$\n",
"x^2(t) + y^2(t) = 1, \\qquad \\forall t\n",
"$$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"import numpy as np\n",
"from matplotlib import pyplot as plt"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 1
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Forward Euler scheme"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def ForwardEuler(h,T):\n",
" N = int(T/h)\n",
" x,y = np.zeros(N),np.zeros(N)\n",
" x[0] = 1.0\n",
" y[0] = 0.0\n",
" for n in range(1,N):\n",
" x[n] = x[n-1] - h*y[n-1]\n",
" y[n] = y[n-1] + h*x[n-1]\n",
" return x,y"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 2
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = ForwardEuler(h,T)\n",
"\n",
"plt.plot(x,y)\n",
"plt.axes().set_aspect('equal')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 3
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The phase space trajectory is spiralling outward."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Backward Euler scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_n = x_{n-1} - h y_n, \\qquad y_n = y_{n-1} + h x_n\n",
"$$\n",
"Eliminate $y_n$ from first equation to get\n",
"$$\n",
"x_n = \\frac{x_{n-1} - h y_{n-1}}{1 + h^2}\n",
"$$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def BackwardEuler(h,T):\n",
" N = int(T/h)\n",
" x,y = np.zeros(N),np.zeros(N)\n",
" x[0] = 1.0\n",
" y[0] = 0.0\n",
" for n in range(1,N):\n",
" x[n] = (x[n-1] - h*y[n-1])/(1.0 + h**2)\n",
" y[n] = y[n-1] + h*x[n]\n",
" return x,y"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 4
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = BackwardEuler(h,T)\n",
"\n",
"plt.plot(x,y)\n",
"plt.axes().set_aspect('equal')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 5
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The phase space trajectory is spiralling inward."
]
},
{
"cell_type": "heading",
"level": 2,
"metadata": {},
"source": [
"Trapezoidal scheme"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"x_n = x_{n-1} - \\frac{h}{2}(y_{n-1} + y_n), \\qquad y_n = y_{n-1} + \\frac{h}{2}(x_{n-1} + x_n)\n",
"$$\n",
"Eliminate $y_n$ from first equation\n",
"$$\n",
"x_n = \\frac{ (1-\\frac{1}{4}h^2) x_{n-1} - h y_{n-1} }{1 + \\frac{1}{4}h^2}\n",
"$$"
]
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"def Trapezoid(h,T):\n",
" N = int(T/h)\n",
" x,y = np.zeros(N),np.zeros(N)\n",
" x[0] = 1.0\n",
" y[0] = 0.0\n",
" for n in range(1,N):\n",
" x[n] = ((1.0-0.25*h**2)*x[n-1] - h*y[n-1])/(1.0 + 0.25*h**2)\n",
" y[n] = y[n-1] + 0.5*h*(x[n-1] + x[n])\n",
" return x,y"
],
"language": "python",
"metadata": {},
"outputs": [],
"prompt_number": 6
},
{
"cell_type": "code",
"collapsed": false,
"input": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = Trapezoid(h,T)\n",
"\n",
"plt.plot(x,y)\n",
"plt.axes().set_aspect('equal')"
],
"language": "python",
"metadata": {},
"outputs": [
{
"metadata": {},
"output_type": "display_data",
"svg": [
"\n",
"\n",
"\n",
"\n"
],
"text": [
""
]
}
],
"prompt_number": 7
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The phase space trajectory is exactly the unit circle. \n",
"\n",
"Multiply first equation by $x_n + x_{n-1}$ and second equation by $y_n + y_{n-1}$\n",
"$$\n",
"(x_n + x_{n-1})(x_n - x_{n-1}) = - \\frac{h}{2}(x_n + x_{n-1})(y_n + y_{n-1})\n",
"$$\n",
"$$\n",
"(y_n + y_{n-1})(y_n - y_{n-1}) = + \\frac{h}{2}(x_n + x_{n-1})(y_n + y_{n-1})\n",
"$$\n",
"Adding the two equations we get\n",
"$$\n",
"x_n^2 + y_n^2 = x_{n-1}^2 + y_{n-1}^2\n",
"$$\n",
"Thus the Trapezoidal method is able to preserve the invariant."
]
}
],
"metadata": {}
}
]
}