{
"cells": [
{
"cell_type": "markdown",
"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",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"import numpy as np\n",
"from matplotlib import pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"theta = np.linspace(0.0, 2*np.pi, 500)\n",
"xe = np.cos(theta)\n",
"ye = np.sin(theta)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forward Euler scheme"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = ForwardEuler(h,T)\n",
"\n",
"plt.plot(xe,ye,'r--',label='Exact')\n",
"plt.plot(x,y,label='FE')\n",
"plt.legend()\n",
"plt.grid(True)\n",
"plt.gca().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The phase space trajectory is spiralling outward."
]
},
{
"cell_type": "markdown",
"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",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = BackwardEuler(h,T)\n",
"\n",
"plt.plot(xe,ye,'r--',label='Exact')\n",
"plt.plot(x,y,label='BE')\n",
"plt.legend()\n",
"plt.grid(True)\n",
"plt.gca().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The phase space trajectory is spiralling inward."
]
},
{
"cell_type": "markdown",
"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",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"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"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = Trapezoid(h,T)\n",
"\n",
"plt.plot(xe,ye,'r--',label='Exact')\n",
"plt.plot(x,y,label='Trap')\n",
"plt.legend()\n",
"plt.grid(True)\n",
"plt.gca().set_aspect('equal')"
]
},
{
"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.\n",
"\n",
"**Excercise**: Show that the energy increases for forward Euler and decreases for backward Euler, for any step size $h$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Partitioned Euler\n",
"This is a symplectic Runge-Kutta method. $x$ is updated explicitly and $y$ is updated implicitly.\n",
"$$\n",
"x_n = x_{n-1} - h y_{n-1}, \\qquad y_n = y_{n-1} + h x_{n}\n",
"$$\n",
"The update of $y$ uses the latest updated value of $x$."
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [],
"source": [
"def PartEuler(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]\n",
" return x,y"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.02\n",
"T = 4.0*np.pi\n",
"x,y = PartEuler(h,T)\n",
"\n",
"plt.plot(xe,ye,'r--',label='Exact')\n",
"plt.plot(x,y,label='PE')\n",
"plt.legend()\n",
"plt.grid(True)\n",
"plt.gca().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We get very good results, even though the method is only first order. We can also switch the implicit/explicit parts\n",
"\n",
"$$\n",
"y_n = y_{n-1} + h x_{n-1}, \\qquad x_n = x_{n-1} - h y_n\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Classical RK4 scheme"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"def rhs(u):\n",
" return np.array([u[1], -u[0]])\n",
"\n",
"def RK4(h,T):\n",
" N = int(T/h)\n",
" u = np.zeros((N,2))\n",
" u[0,0] = 1.0 # x\n",
" u[0,1] = 0.0 # y\n",
" for n in range(0,N-1):\n",
" w = u[n,:]\n",
" k0 = rhs(w)\n",
" k1 = rhs(w + 0.5*h*k0)\n",
" k2 = rhs(w + 0.5*h*k1)\n",
" k3 = rhs(w + h*k2)\n",
" u[n+1,:] = w + (h/6)*(k0 + k1 + k2 + k3)\n",
" return u[:,0], u[:,1]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We test this method for a longer time."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"h = 0.02\n",
"T = 10.0*np.pi\n",
"x,y = RK4(h,T)\n",
"\n",
"plt.plot(xe,ye,'r--',label='Exact')\n",
"plt.plot(x,y,label='RK4')\n",
"plt.legend()\n",
"plt.grid(True)\n",
"plt.gca().set_aspect('equal')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"RK4 is a more accurate method compared to forward/backward Euler schemes, but it still loses total energy."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}