Skip to content

Instantly share code, notes, and snippets.

@victor-cali

victor-cali/numerical_merhods_odes.ipynb

Last active December 12, 2020 01:52
Show Gist options
  • Save victor-cali/34f14e23bfcf321155dd67b3058cf73b to your computer and use it in GitHub Desktop.
Save victor-cali/34f14e23bfcf321155dd67b3058cf73b to your computer and use it in GitHub Desktop.
Download ZIP
Numerical_Merhods_ODEs.ipynb
Display the source blob
Display the rendered blob
Raw
numerical_merhods_odes.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"name": "Numerical_Merhods_ODEs.ipynb",
"provenance": [],
"collapsed_sections": [],
"authorship_tag": "ABX9TyPmtr8tVy2URZFihW/z0CnS",
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/VictorCanLima/34f14e23bfcf321155dd67b3058cf73b/numerical_merhods_odes.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "code",
"metadata": {
"id": "JI5diYcTZx8O"
},
"source": [
"w0 = 0.5\r\n",
"a = 0\r\n",
"b = 2\r\n",
"N = 100\r\n",
"h = (b-a)/N\r\n",
"imax = (b-a)/h\r\n",
"\r\n",
"def f(t,y):\r\n",
" return y-(t**2)+1;"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "ED5USUtPZ2mN"
},
"source": [
"def euler():\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" w_i = w_i+h*f(t_i,w_i)\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def taylor_2nd_order():\r\n",
" def first_derivative(t,w):\r\n",
" return 0\r\n",
" def T(t,w):\r\n",
" return f(t,w)+(h/2)*first_derivative(t,w)\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" w_i = w_i+T(t_i,w_i)\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def taylor_4th_order():\r\n",
" from math import factorial as fact\r\n",
" def first_derivative(t,w):\r\n",
" return 0\r\n",
" def second_derivative(t,w):\r\n",
" return 0\r\n",
" def third_derivative(t,w):\r\n",
" return 0\r\n",
" def T(t,w):\r\n",
" temp1 = f(t,w)+(h/2)*first_derivative(t,w)\r\n",
" temp2 = (pow(h,2)/fact(3))*second_derivative(t,w)\r\n",
" return temp1+temp2+(pow(h,3)/fact(4))*third_derivative(t,w)\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" w_i = w_i+T(t_i,w_i)\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def euler_modified():\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" w_i = w_i+(h/2)*(f(t_i,w_i)+f(t_i+h,w_i+f(t_i,w_i)))\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def runge_kutta_midpoint():\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" w_i = w_i+h*f(t_i+(h/2),w_i+(h/2)*f(t_i,w_i))\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def runge_kutta_order4():\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" for i in range(N):\r\n",
" t_i = t_i + h\r\n",
" k_1 = h*f(t_i,w_i)\r\n",
" k_2 = h*f(t_i+(h/2),w_i+(0.5*k_1))\r\n",
" k_3 = h*f(t_i+(h/2),w_i+(0.5*k_2))\r\n",
" k_4 = h*f(t_i+h,w_i+k_3)\r\n",
" w_i = w_i+(1/6)*(k_1+2*k_2+2*k_3+k_4)\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" return ws, ts\r\n",
"\r\n",
"def runge_kutta_fehldberg(hmin,hmax,tol):\r\n",
" h_fehldberg = hmax\r\n",
" ts = [a]\r\n",
" ws = [w0]\r\n",
" hs = [h_fehldberg]\r\n",
" t_i = a\r\n",
" w_i = w0\r\n",
" flag = True\r\n",
" while flag:\r\n",
" k_1 = h_fehldberg * f(t_i, w_i)\r\n",
" k_2 = h_fehldberg * f(t_i + 0.25*h_fehldberg, w_i + 0.25*k_1)\r\n",
" k_3 = h_fehldberg * f(t_i + 0.375*h_fehldberg, w_i + 0.09375*k_1 + 0.28125*k_2)\r\n",
" k_4 = h_fehldberg * f(t_i + 0.92307692*h_fehldberg, w_i + 0.87938097*k_1 - 3.2771961766*k_2 + 3.32089215*k_3)\r\n",
" k_5 = h_fehldberg * f(t_i + h_fehldberg, w_i + 2.032407407*k_1 - 8*k_2 + 7.17348927*k_3 - 0.205896686*k_4)\r\n",
" k_6 = h_fehldberg * f(t_i + 0.5*h_fehldberg, w_i - 0.29629629*k_1 + 2*k_2 - 1.3816764*k_3 + 0.4533040*k_4 - 0.275*k_5)\r\n",
" R = (1.0/h_fehldberg) * abs(0.0027778*k_1 - 0.0299415 * k_3 - 0.0291998936*k_4 +0.02*k_5 + 0.036363636*k_6)\r\n",
" if R <= tol:\r\n",
" t_i = t_i + h_fehldberg\r\n",
" w_i = w_i + 25.0/216.0*k_1+1408.0/2565.0*k_3+2197.0/4104.0*k_4-0.2*k_5\r\n",
" ts.append(t_i)\r\n",
" ws.append(w_i)\r\n",
" hs.append(h_fehldberg)\r\n",
"\r\n",
" delta = 0.84*pow(float(tol/R),0.25)\r\n",
" if delta <= 0.1:\r\n",
" h_fehldberg = 0.1*h_fehldberg\r\n",
" elif delta >= 4:\r\n",
" h_fehldberg = 4*h_fehldberg\r\n",
" else:\r\n",
" h_fehldberg = delta*h_fehldberg\r\n",
"\r\n",
" if h_fehldberg>hmax:\r\n",
" h_fehldberg = hmax\r\n",
" if t_i >= b:\r\n",
" flag = False\r\n",
" elif t_i+h_fehldberg > b:\r\n",
" h_fehldberg = b-t_i\r\n",
" elif h_fehldberg < hmin:\r\n",
" flag = False\r\n",
" return ws, ts, hs"
],
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"metadata": {
"id": "9UzmtDLEaC6e"
},
"source": [
"#Colors: 'r', 'b', 'g', 'c', 'm', 'y', 'k', 'w'\r\n",
"def plot():\r\n",
" import matplotlib.pyplot as plt\r\n",
" import matplotlib.colors as mcolors\r\n",
" import matplotlib.patches as mpatches\r\n",
" f, ax = plt.subplots()\r\n",
" myhandles = []\r\n",
" #_____________________________________________________________________________\r\n",
" y1,x1 = runge_kutta_order4()\r\n",
" #ax.plot(x1,y1,'r.')\r\n",
" patch1 = mpatches.Patch(color='r', label='RK4')\r\n",
" myhandles.append(patch1)\r\n",
" #_____________________________________________________________________________\r\n",
" y2,x2 = euler()\r\n",
" #ax.plot(x2,y2,'b.')\r\n",
" patch2 = mpatches.Patch(color='b', label='Euler')\r\n",
" myhandles.append(patch2)\r\n",
" #_____________________________________________________________________________\r\n",
" y3,x3 = euler_modified()\r\n",
" ax.plot(x2,y2,'g.')\r\n",
" patch3 = mpatches.Patch(color='g', label='Euler modified')\r\n",
" myhandles.append(patch3)\r\n",
" #_____________________________________________________________________________\r\n",
" y4,x4 = runge_kutta_midpoint()\r\n",
" #ax.plot(x4,y4,'c.')\r\n",
" patch4 = mpatches.Patch(color='c', label='RK midpoint')\r\n",
" myhandles.append(patch4)\r\n",
" #_____________________________________________________________________________\r\n",
" y5,x5,steps = runge_kutta_fehldberg(0.01,0.25,1.0E-5)\r\n",
" ax.plot(x5,y5,'m.')\r\n",
" patch5 = mpatches.Patch(color='m', label='RKF')\r\n",
" myhandles.append(patch5)\r\n",
" #_____________________________________________________________________________\r\n",
" plt.legend(handles=myhandles)\r\n",
" plt.title(\"Solution of the ODE\") \r\n",
" ax.set_xlabel(\"t\")\r\n",
" ax.set_ylabel('w')\r\n",
" plt.show"
],
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment