Created using Colab
This commit is contained in:
udlbook
432
Blogs/BorealisODENumerical.ipynb
Normal file
432
Blogs/BorealisODENumerical.ipynb
Normal file
@@ -0,0 +1,432 @@
|
||||
{
|
||||
"cells": [
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "view-in-github",
|
||||
"colab_type": "text"
|
||||
},
|
||||
"source": [
|
||||
"<a href=\"https://colab.research.google.com/github/udlbook/udlbook/blob/main/Blogs/BorealisODENumerical.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "JXsO7ce7oqeq"
|
||||
},
|
||||
"source": [
|
||||
"# Numerical methods for ODEs\n",
|
||||
"\n",
|
||||
"This blog contains code that accompanies the RBC Borealis blog on numerical methods for ODEs. Contact udlbookmail@gmail.com if you find any problems."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "AnvAKtP_oqes"
|
||||
},
|
||||
"source": [
|
||||
"Import relevant libraries"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "UF-gJyZggyrl"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"import numpy as np\n",
|
||||
"import matplotlib.pyplot as plt"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "szWLVrSSoqet"
|
||||
},
|
||||
"source": [
|
||||
"Define the ODE that we will be experimenting with."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "NkrGZLL6iM3P"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# The ODE that we will experiment with\n",
|
||||
"def ode_lin_homog(t,x):\n",
|
||||
" return 0.5 * x ;\n",
|
||||
"\n",
|
||||
"# The derivative of the ODE function with respect to x (needed for Taylor's method)\n",
|
||||
"def ode_lin_homog_deriv_x(t,x):\n",
|
||||
" return 0.5 ;\n",
|
||||
"\n",
|
||||
"# The derivative of the ODE function with respect to t (needed for Taylor's method)\n",
|
||||
"def ode_lin_homog_deriv_t(t,x):\n",
|
||||
" return 0.0 ;\n",
|
||||
"\n",
|
||||
"# The closed form solution (so we can measure the error)\n",
|
||||
"def ode_lin_homog_soln(t,C=0.5):\n",
|
||||
" return C * np.exp(0.5 * t) ;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "In1C9wZkoqet"
|
||||
},
|
||||
"source": [
|
||||
"This is a generic method that runs the numerical methods. It takes the initial conditions ($t_0$, $x_0$), the final time $t_1$ and the step size $h$. It also takes the ODE function itself and its derivatives (only used for Taylor's method). Finally, the parameter \"step_function\" is the method used to update (e.g., Euler's methods, Runge-Kutte 4-step)."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "VZfZDJAfmyrf"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def run_numerical(x_0, t_0, t_1, h, ode_func, ode_func_deriv_x, ode_func_deriv_t, ode_soln, step_function):\n",
|
||||
" x = [x_0]\n",
|
||||
" t = [t_0]\n",
|
||||
" while (t[-1] <= t_1):\n",
|
||||
" x = x+[step_function(x[-1],t[-1],h, ode_func, ode_func_deriv_x, ode_func_deriv_t)]\n",
|
||||
" t = t + [t[-1]+h]\n",
|
||||
"\n",
|
||||
" # Returns x,y plot plus total numerical error at last point.\n",
|
||||
" return t, x, np.abs(ode_soln(t[-1])-x[-1])"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "Vfkc3-_7oqet"
|
||||
},
|
||||
"source": [
|
||||
"Run the numerical method with step sizes of 2.0, 1.0, 0.5, 0.25, 0.125, 0.0675 and plot the results"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "1tyGbMZhoqeu"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def run_and_plot(ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function):\n",
|
||||
" # Specify the grid of points to draw the ODE\n",
|
||||
" t = np.arange(0.04, 4.0, 0.2)\n",
|
||||
" x = np.arange(0.04, 4.0, 0.2)\n",
|
||||
" T, X = np.meshgrid(t,x)\n",
|
||||
"\n",
|
||||
" # ODE equation at these grid points (used to draw quiver-plot)\n",
|
||||
" dx = ode(T,X)\n",
|
||||
" dt = np.ones(dx.shape)\n",
|
||||
"\n",
|
||||
" # The ground truth solution\n",
|
||||
" t2= np.arange(0,10,0.1)\n",
|
||||
" x2 = ode_solution(t2)\n",
|
||||
"\n",
|
||||
" #####################################x_0, t_0, t_1, h #################################################\n",
|
||||
" t_sim1,x_sim1,error1 = run_numerical(0.5, 0.0, 4.0, 2.0000, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" t_sim2,x_sim2,error2 = run_numerical(0.5, 0.0, 4.0, 1.0000, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" t_sim3,x_sim3,error3 = run_numerical(0.5, 0.0, 4.0, 0.5000, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" t_sim4,x_sim4,error4 = run_numerical(0.5, 0.0, 4.0, 0.2500, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" t_sim5,x_sim5,error5 = run_numerical(0.5, 0.0, 4.0, 0.1250, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" t_sim6,x_sim6,error6 = run_numerical(0.5, 0.0, 4.0, 0.0675, ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
"\n",
|
||||
" # Plot the ODE and ground truth solution\n",
|
||||
" fig,ax = plt.subplots()\n",
|
||||
" ax.quiver(T,X,dt,dx, scale=35.0)\n",
|
||||
" ax.plot(t2,x2,'r-')\n",
|
||||
"\n",
|
||||
" # Plot the numerical approximations\n",
|
||||
" ax.plot(t_sim1,x_sim1,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
" ax.plot(t_sim2,x_sim2,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
" ax.plot(t_sim3,x_sim3,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
" ax.plot(t_sim4,x_sim4,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
" ax.plot(t_sim5,x_sim5,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
" ax.plot(t_sim6,x_sim6,'.-',markeredgecolor='#773c23ff',markerfacecolor='#d18362', color='#d18362', markersize=10)\n",
|
||||
"\n",
|
||||
" ax.set_aspect('equal')\n",
|
||||
" ax.set_xlim(0,4)\n",
|
||||
" ax.set_ylim(0,4)\n",
|
||||
"\n",
|
||||
" plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "JYrq8QIwvOIy"
|
||||
},
|
||||
"source": [
|
||||
"# Euler Method\n",
|
||||
"\n",
|
||||
"Define the Euler method and set up functions for plotting."
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "N73xMnCukVVX"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def euler_step(x_0, t_0, h, ode_func, ode_func_deriv_x=None, ode_func_deriv_t=None):\n",
|
||||
" return x_0 + h * ode_func(t_0, x_0) ;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "4B1_PGEcsZ9H"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"run_and_plot(ode_lin_homog, None, None, ode_lin_homog_soln, euler_step)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "FfwNihtkvJeX"
|
||||
},
|
||||
"source": [
|
||||
"# Heun's Method"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "srHfNDcDxI1o"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def heun_step(x_0, t_0, h, ode_func, ode_func_deriv_x=None, ode_func_deriv_t=None):\n",
|
||||
" f_x0_t0 = ode_func(t_0, x_0)\n",
|
||||
" return x_0 + h/2 * ( f_x0_t0 + ode_func(t_0+h, x_0+h*f_x0_t0)) ;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "WOApHz9xoqev"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"run_and_plot(ode_lin_homog, None, None, ode_lin_homog_soln, heun_step)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "0XSzzFDIvRhm"
|
||||
},
|
||||
"source": [
|
||||
"# Modified Euler method"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "fSXprgVJ5Yep"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def modified_euler_step(x_0, t_0, h, ode_func, ode_func_deriv_x=None, ode_func_deriv_t=None):\n",
|
||||
" f_x0_t0 = ode_func(t_0, x_0)\n",
|
||||
" return x_0 + h * ode_func(t_0+h/2, x_0+ h * f_x0_t0/2) ;"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "8LKSrCD2oqev"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"run_and_plot(ode_lin_homog, None, None, ode_lin_homog_soln, modified_euler_step)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "yp8ZBpwooqev"
|
||||
},
|
||||
"source": [
|
||||
"# Second order Taylor's method"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "NtBBgzWLoqev"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def taylor_2nd_order(x_0, t_0, h, ode_func, ode_func_deriv_x, ode_func_deriv_t):\n",
|
||||
" f1 = ode_func(t_0, x_0)\n",
|
||||
" return x_0 + h * ode_func(t_0, x_0) + (h*h/2) * (ode_func_deriv_x(t_0,x_0) * ode_func(t_0, x_0) + ode_func_deriv_t(t_0, x_0))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "ioeeIohUoqev"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"run_and_plot(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, taylor_2nd_order)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "WcuhV5lL1zAJ"
|
||||
},
|
||||
"source": [
|
||||
"# Fourth Order Runge Kutta"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "0NZN81Bpwu56"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def runge_kutta_4_step(x_0, t_0, h, ode_func, ode_func_deriv_x=None, ode_func_deriv_t=None):\n",
|
||||
" f1 = ode_func(t_0, x_0)\n",
|
||||
" f2 = ode_func(t_0+h/2,x_0+f1 * h/2)\n",
|
||||
" f3 = ode_func(t_0+h/2,x_0+f2 * h/2)\n",
|
||||
" f4 = ode_func(t_0+h, x_0+ f3*h)\n",
|
||||
" return x_0 + (h/6) * (f1 + 2*f2 + 2*f3+f4)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "K-OxE9E6oqew"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"run_and_plot(ode_lin_homog, None, None, ode_lin_homog_soln, runge_kutta_4_step)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "7JifxBhhoqew"
|
||||
},
|
||||
"source": [
|
||||
"# Plot the error as a function of step size"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "ZoEpmlCfsi9P"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Run systematically with a number of different step sizes and store errors for each\n",
|
||||
"def get_errors(ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function):\n",
|
||||
" # Choose the step size h to divide the plotting interval into 1,2,4,8... segments.\n",
|
||||
" # The plots in the article add a few more smaller step sizes, but this takes a while to compute.\n",
|
||||
" # Add them back in if you want the full plot.\n",
|
||||
" all_h = (1./np.array([1,2,4,8,16,32,64,128,256,512,1024,2048,4096])).tolist()\n",
|
||||
" all_err = []\n",
|
||||
"\n",
|
||||
" for i in range(len(all_h)):\n",
|
||||
" t_sim,x_sim,err = run_numerical(0.5, 0.0, 4.0, all_h[i], ode, ode_deriv_x, ode_deriv_t, ode_solution, step_function)\n",
|
||||
" all_err = all_err + [err]\n",
|
||||
"\n",
|
||||
" return all_h, all_err"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {
|
||||
"id": "X0O0KK47xF28"
|
||||
},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# Plot the errors\n",
|
||||
"all_h, all_err_euler = get_errors(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, euler_step)\n",
|
||||
"all_h, all_err_heun = get_errors(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, heun_step)\n",
|
||||
"all_h, all_err_mod_euler = get_errors(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, modified_euler_step)\n",
|
||||
"all_h, all_err_taylor = get_errors(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, taylor_2nd_order)\n",
|
||||
"all_h, all_err_rk = get_errors(ode_lin_homog, ode_lin_homog_deriv_x, ode_lin_homog_deriv_t, ode_lin_homog_soln, runge_kutta_4_step)\n",
|
||||
"\n",
|
||||
"\n",
|
||||
"fig, ax = plt.subplots()\n",
|
||||
"ax.loglog(all_h, all_err_euler,'ro-')\n",
|
||||
"ax.loglog(all_h, all_err_heun,'bo-')\n",
|
||||
"ax.loglog(all_h, all_err_mod_euler,'go-')\n",
|
||||
"ax.loglog(all_h, all_err_taylor,'co-')\n",
|
||||
"ax.loglog(all_h, all_err_rk,'mo-')\n",
|
||||
"ax.set_ylim(1e-13,1e1)\n",
|
||||
"ax.set_xlim(1e-6,1e1)\n",
|
||||
"ax.set_aspect(0.5)\n",
|
||||
"ax.set_xlabel('Step size, $h$')\n",
|
||||
"ax.set_ylabel('Error')\n",
|
||||
"plt.show()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {
|
||||
"id": "BttOqpeo9MsJ"
|
||||
},
|
||||
"source": [
|
||||
"Note that for this ODE, the Heun, Modified Euler and Taylor methods provide EXACTLY the same updates, and so the error curves for all three are identical (subject to difference is numerical rounding errors). This is not in general the case, although the general trend would be the same for each."
|
||||
]
|
||||
}
|
||||
],
|
||||
"metadata": {
|
||||
"colab": {
|
||||
"provenance": [],
|
||||
"include_colab_link": true
|
||||
},
|
||||
"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.9.10"
|
||||
}
|
||||
},
|
||||
"nbformat": 4,
|
||||
"nbformat_minor": 0
|
||||
}
|
||||
Reference in New Issue
Block a user