diff --git a/18_initial_value_ode/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb b/18_initial_value_ode/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb
deleted file mode 100644
index 2306605..0000000
--- a/18_initial_value_ode/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb
+++ /dev/null
@@ -1,837 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {
- "collapsed": true,
- "slideshow": {
- "slide_type": "skip"
- }
- },
- "outputs": [],
- "source": [
- "%plot --format svg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {
- "collapsed": true,
- "slideshow": {
- "slide_type": "skip"
- }
- },
- "outputs": [],
- "source": [
- "setdefaults"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "# Initial Value Problems (ODEs)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "# Initial value problems are integrals\n",
- "\n",
- "$y'+2y=0$\n",
- "\n",
- "$y(0)=1$\n",
- "\n",
- "\n",
- "Solve for $y(t)$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "fragment"
- }
- },
- "source": [
- "$\\frac{dy}{dt}=-2y$\n",
- "\n",
- "$\\frac{dy}{y}=-2dt$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "fragment"
- }
- },
- "source": [
- "Integrate both sides:\n",
- "\n",
- "$\\ln{\\frac{y}{y_{0}}}=-2t+2t_{0}$\n",
- "\n",
- "$y(t)=y_{0}e^{-2t+2t_{0}}=e^{-2t}$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "subslide"
- }
- },
- "source": [
- "## Euler's method \n",
- "\n",
- "$\\frac{dy}{dt}=f(t,y)$\n",
- "\n",
- "$y_{i+1}=y_{i}+\\int_{t_{i}}^{t_{i+1}}f(t,y)dt$\n",
- "\n",
- "$y_{i+1}\\approx y_{i}+f(t_{i},y_{i})h$\n",
- "\n",
- "The error of this method is:\n",
- "\n",
- "$E_{t}=\\frac{f'(t_i , y_i )}{2!}h^2 + \\cdots + O(h^{n+1})$\n",
- "\n",
- "or\n",
- "\n",
- "$E_{a}=O(h^2)$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "### Example: Freefalling problem\n",
- "\n",
- "An object is falling and has a drag coefficient of 0.25 kg/m and mass of 60 kg\n",
- "Define time from 0 to 12 seconds with `N` timesteps \n",
- "function defined as `freefall`\n",
- "\n",
- "Using the Euler ODE solution results in a conditionally stable solution *(at some point the time steps are too large to solve the problem)*"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "subslide"
- }
- },
- "outputs": [],
- "source": [
- "function [v_analytical,v_terminal,t]=freefall(N,tmax)\n",
- " t=linspace(0,tmax,N)';\n",
- " c=0.25; m=60; g=9.81; v_terminal=sqrt(m*g/c);\n",
- "\n",
- " v_analytical = v_terminal*tanh(g*t/v_terminal);\n",
- " v_numerical=zeros(length(t),1);\n",
- " delta_time =diff(t);\n",
- " for i=1:length(t)-1\n",
- " v_numerical(i+1)=v_numerical(i)+(g-c/m*v_numerical(i)^2)*delta_time(i);\n",
- " end\n",
- " % Print values near 0,2,4,6,8,10,12 seconds\n",
- " indices = round(linspace(1,length(t),7));\n",
- " fprintf('time (s)| error (m/s)\\n')\n",
- " fprintf('-------------------------\\n')\n",
- " M=[t(indices),abs(v_analytical(indices)-v_numerical(indices))];\n",
- " fprintf('%7.1f | %10.2f\\n',M(:,1:2)');\n",
- " plot(t,v_analytical,'-',t,v_numerical,'o-')\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "subslide"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "time (s)| error (m/s)\n",
- "-------------------------\n",
- " 0.0 | 0.00\n",
- " 2.3 | 0.46\n",
- " 4.0 | 0.95\n",
- " 6.3 | 1.03\n",
- " 8.0 | 0.80\n",
- " 10.3 | 0.46\n",
- " 12.0 | 0.28\n",
- "\n",
- "O(h^2)=0.33\n"
- ]
- },
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "[v_an,v_t,t]=freefall(22,12);\n",
- "fprintf('\\nO(h^2)=%1.2f',min(diff(t).^2))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "## Heun Method\n",
- "\n",
- "Increase accuracy with *predictor-corrector approach*\n",
- "\n",
- "$y_{i+1}=y_{i}^{m}+f(t_{i},y_{i})h$\n",
- "\n",
- "$y_{i+1}^{j}=y_{i}^{m}+\n",
- "\\frac{f(t_{i},y_{i}^{m})+f(t_{i+1},y_{i+1}^{i-1})}{2}h$\n",
- "\n",
- "This is analagous to the trapezoidal rule\n",
- "\n",
- "$\\int_{t_{i}}^{t_{i+1}}f(t,y)dt=\\frac{f(t_{i},y_{i})+f(t_{i+1},y_{i+1})}{2}h$\n",
- "\n",
- "therefore the error is\n",
- "\n",
- "$E_{t}=\\frac{-f''(\\xi)}{12}h^3$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "subslide"
- }
- },
- "source": [
- "### Example with Heun's method\n",
- "\n",
- "Problem Statement. Use Heun’s method with iteration to integrate \n",
- "\n",
- "$y' = 4e^{0.8t} − 0.5y$\n",
- "\n",
- "from t = 0 to 4 with a step size of 1. The initial condition at t = 0 is y = 2. Employ a stopping criterion of 0.00001% to terminate the corrector iterations."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 3,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "subslide"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "dy =\n",
- "\n",
- " 3\n",
- " 0\n",
- " 0\n",
- " 0\n",
- " 0\n",
- "\n",
- "y =\n",
- "\n",
- " 2\n",
- " 5\n",
- " 0\n",
- " 0\n",
- " 0\n",
- "\n"
- ]
- }
- ],
- "source": [
- "yp=@(t,y) 4*exp(0.8*t)-0.5*y;\n",
- "t=linspace(0,4,5)';\n",
- "y=zeros(size(t));\n",
- "dy=zeros(size(t));\n",
- "dy_corr=zeros(size(t));\n",
- "y(1)=2;\n",
- "dy(1)=yp(t(1),y(1))\n",
- "y(2)=y(1)+dy(1)*(t(2)-t(1))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 4,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "dy_corr =\n",
- "\n",
- " 4.70108\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n",
- "y =\n",
- "\n",
- " 2.00000\n",
- " 6.70108\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n"
- ]
- }
- ],
- "source": [
- "% improve estimate for y(2)\n",
- "dy_corr(1)=(dy(1)+yp(t(2),y(2)))/2\n",
- "y(2)=y(1)+dy_corr(1)*(t(2)-t(1))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "### This process can be iterated until a desired tolerance is achieved"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "fragment"
- }
- },
- "outputs": [],
- "source": [
- "yp=@(t,y) 4*exp(0.8*t)-0.5*y;\n",
- "t=linspace(0,4,5)';\n",
- "y=zeros(size(t));\n",
- "dy=zeros(size(t));\n",
- "dy_corr=zeros(size(t));\n",
- "y(1)=2;\n",
- "for i=1:length(t)-1\n",
- " dy(i)=yp(t(i),y(i));\n",
- " dy_corr(i)=yp(t(i),y(i));\n",
- " y(i+1)=y(i)+dy_corr(i)*(t(i+1)-t(i));\n",
- " n=0;\n",
- " e=10;\n",
- " while (1)\n",
- " n=n+1;\n",
- " yold=y(i+1);\n",
- " dy_corr(i)=(dy(i)+yp(t(i+1),y(i+1)))/2;\n",
- " y(i+1)=y(i)+dy_corr(i)*(t(i+1)-t(i));\n",
- " e=abs(y(i+1)-yold)/y(i+1)*100;\n",
- " if e<= 0.00001 | n>100, break, end\n",
- " end\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 14,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "y_an =\n",
- "\n",
- "@(t) 4 / 1.3 * exp (0.8 * t) - 1.0769 * exp (-t / 2)\n",
- "\n",
- "dy_an =\n",
- "\n",
- "@(t) 0.8 * 4 / 1.3 * exp (0.8 * t) + 1.0769 / 2 * exp (-t / 2)\n",
- "\n"
- ]
- }
- ],
- "source": [
- "\n",
- "y_euler=zeros(size(t));\n",
- "for i=1:length(t)-1\n",
- " dy(i)=yp(t(i),y(i));\n",
- " y_euler(i+1)=y_euler(i)+dy(i)*(t(i+1)-t(i));\n",
- "end\n",
- "\n",
- "y_an =@(t) 4/1.3*exp(0.8*t)-1.0769*exp(-t/2)\n",
- "dy_an=@(t) 0.8*4/1.3*exp(0.8*t)+1.0769/2*exp(-t/2)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 15,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "plot(t,y,'o',t,y_euler,'s',linspace(min(t),max(t)),y_an(linspace(min(t),max(t))))\n",
- "legend('Heuns method','Euler','analytical','Location','NorthWest')\n",
- "xlabel('time')\n",
- "ylabel('y')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": null,
- "metadata": {
- "collapsed": true
- },
- "outputs": [],
- "source": []
- }
- ],
- "metadata": {
- "celltoolbar": "Slideshow",
- "kernelspec": {
- "display_name": "Octave",
- "language": "octave",
- "name": "octave"
- },
- "language_info": {
- "file_extension": ".m",
- "help_links": [
- {
- "text": "MetaKernel Magics",
- "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
- }
- ],
- "mimetype": "text/x-octave",
- "name": "octave",
- "version": "0.19.14"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/18_initial_value_ode/.ipynb_checkpoints/lecture_22-checkpoint.ipynb b/18_initial_value_ode/.ipynb_checkpoints/lecture_22-checkpoint.ipynb
deleted file mode 100644
index 8fdd6a6..0000000
--- a/18_initial_value_ode/.ipynb_checkpoints/lecture_22-checkpoint.ipynb
+++ /dev/null
@@ -1,821 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 1,
- "metadata": {
- "collapsed": true
- },
- "outputs": [],
- "source": [
- "%plot --format svg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 2,
- "metadata": {
- "collapsed": true
- },
- "outputs": [],
- "source": [
- "setdefaults"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n",
- "\n"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "# Initial Value Problem\n",
- "\n",
- "## Euler's method\n",
- "\n",
- "$\\frac{dy}{dt}=f(t,y)$\n",
- "\n",
- "$y_{i+1}=y_{i}+f(t_{i},y_{i})h$\n",
- "\n",
- "The error of this method is:\n",
- "\n",
- "$E_{t}=\\frac{f'(t_i , y_i )}{2!}h^2 + \\cdots + O(h^{n+1})$\n",
- "\n",
- "or\n",
- "\n",
- "$E_{a}=O(h^2)$\n",
- "\n",
- "### Example: Freefalling problem\n",
- "\n",
- "An object is falling and has a drag coefficient of 0.25 kg/m and mass of 60 kg\n",
- "Define time from 0 to 12 seconds with `N` timesteps \n",
- "function defined as `freefall`\n",
- "\n",
- "Using the Euler ODE solution results in a conditionally stable solution *(at some point the time steps are too large to solve the problem)*"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "function [v_analytical,v_terminal,t]=freefall(N,tmax)\n",
- " t=linspace(0,tmax,N)';\n",
- " c=0.25; m=60; g=9.81; v_terminal=sqrt(m*g/c);\n",
- "\n",
- " v_analytical = v_terminal*tanh(g*t/v_terminal);\n",
- " v_numerical=zeros(length(t),1);\n",
- " delta_time =diff(t);\n",
- " for i=1:length(t)-1\n",
- " v_numerical(i+1)=v_numerical(i)+(g-c/m*v_numerical(i)^2)*delta_time(i);\n",
- " end\n",
- " % Print values near 0,2,4,6,8,10,12 seconds\n",
- " indices = round(linspace(1,length(t),7));\n",
- " fprintf('time (s)| error (m/s)\\n')\n",
- " fprintf('-------------------------\\n')\n",
- " M=[t(indices),abs(v_analytical(indices)-v_numerical(indices))];\n",
- " fprintf('%7.1f | %10.2f\\n',M(:,1:2)');\n",
- " plot(t,v_analytical,'-',t,v_numerical,'o-')\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 8,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "time (s)| error (m/s)\n",
- "-------------------------\n",
- " 0.0 | 0.00\n",
- " 7.1 | 26.67\n",
- " 14.3 | 54.21\n",
- " 28.6 | 33.62\n",
- " 35.7 | 29.84\n",
- " 42.9 | 82.85\n",
- " 50.0 | 47.86\n",
- "\n",
- "O(h^2)=51.02\n"
- ]
- },
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "[v_an,v_t,t]=freefall(8,50);\n",
- "fprintf('\\nO(h^2)=%1.2f',min(diff(t).^2))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "## Heun Method\n",
- "\n",
- "Increase accuracy with *predictor-corrector approach*\n",
- "\n",
- "$y_{i+1}=y_{i}^{m}+f(t_{i},y_{i})h$\n",
- "\n",
- "$y_{i+1}^{j}=y_{i}^{m}+\n",
- "\\frac{f(t_{i},y_{i}^{m})+f(t_{i+1},y_{i+1}^{i-1})}{2}h$\n",
- "\n",
- "This is analagous to the trapezoidal rule\n",
- "\n",
- "$\\int_{t_{i}}^{t_{i+1}}f(t)dt=\\frac{f(t_{i})+f(t_{i+1})}{2}h$\n",
- "\n",
- "therefore the error is\n",
- "\n",
- "$E_{t}=\\frac{-f''(\\xi)}{12}h^3$"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "**Example with Heun's method**\n",
- "\n",
- "Problem Statement. Use Heun’s method with iteration to integrate \n",
- "\n",
- "$y' = 4e^{0.8t} − 0.5y$\n",
- "\n",
- "from t = 0 to 4 with a step size of 1. The initial condition at t = 0 is y = 2. Employ a stopping criterion of 0.00001% to terminate the corrector iterations."
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 102,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "dy =\n",
- "\n",
- " 3\n",
- " 0\n",
- " 0\n",
- " 0\n",
- " 0\n",
- "\n",
- "y =\n",
- "\n",
- " 2\n",
- " 5\n",
- " 0\n",
- " 0\n",
- " 0\n",
- "\n"
- ]
- }
- ],
- "source": [
- "yp=@(t,y) 4*exp(0.8*t)-0.5*y;\n",
- "t=linspace(0,4,5)';\n",
- "y=zeros(size(t));\n",
- "dy=zeros(size(t));\n",
- "dy_corr=zeros(size(t));\n",
- "y(1)=2;\n",
- "dy(1)=yp(t(1),y(1))\n",
- "y(2)=y(1)+dy(1)*(t(2)-t(1))"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 103,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "dy_corr =\n",
- "\n",
- " 4.70108\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n",
- "y =\n",
- "\n",
- " 2.00000\n",
- " 6.70108\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n"
- ]
- }
- ],
- "source": [
- "% improve estimate for y(2)\n",
- "dy_corr(1)=(dy(1)+yp(t(2),y(2)))/2\n",
- "y(2)=y(1)+dy_corr(1)*(t(2)-t(1))"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {},
- "source": [
- "This process can be iterated until a desired tolerance is achieved"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 106,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "n=0;\n",
- "e=10;\n",
- "while (1)\n",
- " n=n+1;\n",
- " yold=y(2);\n",
- " dy_corr(1)=(dy(1)+yp(t(2),y(2)))/2;\n",
- " y(2)=y(1)+dy_corr(1)*(t(2)-t(1));\n",
- " e=abs(y(2)-yold)/y(2)*100;\n",
- " if e<= 0.00001 | n>100, break, end\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 107,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "y =\n",
- "\n",
- " 2.00000\n",
- " 6.36087\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n",
- "dy_corr =\n",
- "\n",
- " 4.36087\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- " 0.00000\n",
- "\n",
- "n = 12\n",
- "e = 6.3760e-06\n"
- ]
- }
- ],
- "source": [
- "y\n",
- "dy_corr\n",
- "n\n",
- "e"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 127,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "y_an =\n",
- "\n",
- "@(t) 4 / 1.3 * exp (0.8 * t) - 1.0769 * exp (-t / 2)\n",
- "\n",
- "dy_an =\n",
- "\n",
- "@(t) 0.8 * 4 / 1.3 * exp (0.8 * t) + 1.0769 / 2 * exp (-t / 2)\n",
- "\n"
- ]
- }
- ],
- "source": [
- "\n",
- "y_euler=zeros(size(t));\n",
- "for i=1:length(t)-1\n",
- " y_euler(i+1)=y_euler(i)+dy(i)*(t(i+1)-t(i));\n",
- "end\n",
- "\n",
- "y_an =@(t) 4/1.3*exp(0.8*t)-1.0769*exp(-t/2)\n",
- "dy_an=@(t) 0.8*4/1.3*exp(0.8*t)+1.0769/2*exp(-t/2)\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 121,
- "metadata": {
- "collapsed": false
- },
- "outputs": [],
- "source": [
- "yp=@(t,y) 4*exp(0.8*t)-0.5*y;\n",
- "t=linspace(0,4,5)';\n",
- "y=zeros(size(t));\n",
- "dy=zeros(size(t));\n",
- "dy_corr=zeros(size(t));\n",
- "y(1)=2;\n",
- "for i=1:length(t)-1\n",
- " dy(i)=yp(t(i),y(i));\n",
- " dy_corr(i)=yp(t(i),y(i));\n",
- " y(i+1)=y(i)+dy_corr(i)*(t(i+1)-t(i));\n",
- " n=0;\n",
- " e=10;\n",
- " while (1)\n",
- " n=n+1;\n",
- " yold=y(i+1);\n",
- " dy_corr(i)=(dy(i)+yp(t(i+1),y(i+1)))/2;\n",
- " y(i+1)=y(i)+dy_corr(i)*(t(i+1)-t(i));\n",
- " e=abs(y(i+1)-yold)/y(i+1)*100;\n",
- " if e<= 0.00001 | n>100, break, end\n",
- " end\n",
- "end"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 122,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "y =\n",
- "\n",
- " 2.0000\n",
- " 6.3609\n",
- " 15.3022\n",
- " 34.7433\n",
- " 77.7351\n",
- "\n",
- "ans =\n",
- "\n",
- " 2.0000\n",
- " 6.1946\n",
- " 14.8439\n",
- " 33.6772\n",
- " 75.3390\n",
- "\n"
- ]
- }
- ],
- "source": [
- "y\n",
- "y_an(t)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 128,
- "metadata": {
- "collapsed": false
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "plot(t,y,'o',t,y_euler,'s',linspace(min(t),max(t)),y_an(linspace(min(t),max(t))))\n",
- "legend('Heuns method','Euler','analytical','Location','NorthWest')\n",
- "xlabel('time')\n",
- "ylabel('y')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "collapsed": true
- },
- "source": [
- "# Phugoid Oscillation example\n",
- "\n",
- "[phugoid lessons in python]("
- ]
- }
- ],
- "metadata": {
- "kernelspec": {
- "display_name": "Octave",
- "language": "octave",
- "name": "octave"
- },
- "language_info": {
- "file_extension": ".m",
- "help_links": [
- {
- "text": "MetaKernel Magics",
- "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md"
- }
- ],
- "mimetype": "text/x-octave",
- "name": "octave",
- "version": "0.19.14"
- }
- },
- "nbformat": 4,
- "nbformat_minor": 2
-}
diff --git a/18_initial_value_ode/octave-workspace b/18_initial_value_ode/octave-workspace
deleted file mode 100644
index 8c437bb..0000000
Binary files a/18_initial_value_ode/octave-workspace and /dev/null differ
diff --git a/19_nonlinear_ivp/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb b/19_nonlinear_ivp/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb
deleted file mode 100644
index 361a7bc..0000000
--- a/19_nonlinear_ivp/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb
+++ /dev/null
@@ -1,3482 +0,0 @@
-{
- "cells": [
- {
- "cell_type": "code",
- "execution_count": 11,
- "metadata": {
- "collapsed": true,
- "slideshow": {
- "slide_type": "skip"
- }
- },
- "outputs": [],
- "source": [
- "setdefaults"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 12,
- "metadata": {
- "collapsed": true,
- "slideshow": {
- "slide_type": "skip"
- }
- },
- "outputs": [],
- "source": [
- "%plot --format svg"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 13,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "skip"
- }
- },
- "outputs": [],
- "source": [
- "pkg load odepkg"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "# Initial Value Problems (continued)\n",
- "*ch. 23 - Adaptive methods*"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "## Predator-Prey Models (Chaos Theory)"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "fragment"
- }
- },
- "source": [
- "Predator-prey models were developed independently in the early part of the twentieth\n",
- "century by the Italian mathematician Vito Volterra and the American biologist Alfred\n",
- "Lotka. These equations are commonly called Lotka-Volterra equations. The simplest version is the following pairs of ODEs:\n",
- "\n",
- "$\\frac{dx}{dt}=ax-bxy$\n",
- "\n",
- "$\\frac{dy}{dt}=-cy+dxy$\n",
- "\n",
- "where x and y = the number of prey and predators, respectively, a = the prey growth rate, c = the predator death rate, and b and d = the rates characterizing the effect of the predator-prey interactions on the prey death and the predator growth, respectively."
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "```matlab\n",
- "function yp=predprey(t,y,a,b,c,d)\n",
- " % predator-prey model (Lotka-Volterra equations)\n",
- " yp=zeros(1,2);\n",
- " x=y(1); % population in thousands of prey\n",
- " y=y(2); % population in thousands of predators\n",
- " yp(1)=a.*x-b.*x.*y; % population change in thousands of prey/year\n",
- " yp(2)=-c*y+d*x.*y; % population change in thousands of predators/year\n",
- "end\n",
- "```"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 14,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "fragment"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "ans =\n",
- "\n",
- " 0\n",
- " 19\n",
- "\n"
- ]
- }
- ],
- "source": [
- "predprey(0,[20 1],1,1,1,1)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 15,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "h=0.000625;tspan=[0 40];y0=[2 1];\n",
- "a=1.2;b=0.6;c=0.8;d=0.3;\n",
- "[t y] = eulode(@predprey,tspan,y0,h,a,b,c,d);\n",
- "\n",
- "plot(t,y(:,1),t,y(:,2),'--')\n",
- "legend('prey','predator');\n",
- "title('Euler time plot')\n",
- "xlabel('time (years)')\n",
- "ylabel('population in thousands')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 16,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "plot(y(:,1),y(:,2))\n",
- "title('Euler phase plane plot')\n",
- "xlabel('thousands of prey')\n",
- "ylabel('thousands of predators')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "# Fourth Order Runge-Kutte integration\n",
- "\n",
- "```matlab\n",
- "function [tp,yp] = rk4sys(dydt,tspan,y0,h,varargin)\n",
- "% rk4sys: fourth-order Runge-Kutta for a system of ODEs\n",
- "% [t,y] = rk4sys(dydt,tspan,y0,h,p1,p2,...): integrates\n",
- "% a system of ODEs with fourth-order RK method\n",
- "% input:\n",
- "% dydt = name of the M-file that evaluates the ODEs\n",
- "% tspan = [ti, tf]; initial and final times with output\n",
- "% generated at interval of h, or\n",
- "% = [t0 t1 ... tf]; specific times where solution output\n",
- "% y0 = initial values of dependent variables\n",
- "% h = step size\n",
- "% p1,p2,... = additional parameters used by dydt\n",
- "% output:\n",
- "% tp = vector of independent variable\n",
- "% yp = vector of solution for dependent variables\n",
- "if nargin<4,error('at least 4 input arguments required'), end\n",
- "if any(diff(tspan)<=0),error('tspan not ascending order'), end\n",
- "n = length(tspan);\n",
- "ti = tspan(1);tf = tspan(n);\n",
- "if n == 2\n",
- " t = (ti:h:tf)'; n = length(t);\n",
- " if t(n)h,hh = h;end\n",
- " while(1)\n",
- " if tt+hh>tend,hh = tend-tt;end\n",
- " k1 = dydt(tt,y(i,:),varargin{:})';\n",
- " ymid = y(i,:) + k1.*hh./2;\n",
- " k2 = dydt(tt+hh/2,ymid,varargin{:})';\n",
- " ymid = y(i,:) + k2*hh/2;\n",
- " k3 = dydt(tt+hh/2,ymid,varargin{:})';\n",
- " yend = y(i,:) + k3*hh;\n",
- " k4 = dydt(tt+hh,yend,varargin{:})';\n",
- " phi = (k1+2*(k2+k3)+k4)/6;\n",
- " y(i+1,:) = y(i,:) + phi*hh;\n",
- " tt = tt+hh;\n",
- " i=i+1;\n",
- " if tt>=tend,break,end\n",
- " end\n",
- " np = np+1; tp(np) = tt; yp(np,:) = y(i,:);\n",
- " if tt>=tf,break,end\n",
- "end\n",
- "```"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 17,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "h=0.0625;tspan=[0 40];y0=[2 1];\n",
- "a=1.2;b=0.6;c=0.8;d=0.3;\n",
- "tspan=[0 10];\n",
- "[t y] = rk4sys(@predprey,tspan,y0,h,a,b,c,d);\n",
- "plot(t,y(:,1),t,y(:,2),'--')\n",
- "title('RK4 time plot')\n",
- "xlabel('time (years)')\n",
- "ylabel('population in thousands')"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 18,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- ""
- ],
- "text/plain": [
- ""
- ]
- },
- "metadata": {},
- "output_type": "display_data"
- }
- ],
- "source": [
- "plot(y(:,1),y(:,2))\n",
- "title('RK4 phase plane plot')\n",
- "xlabel('thousands of prey')\n",
- "ylabel('thousands of predators')"
- ]
- },
- {
- "cell_type": "markdown",
- "metadata": {
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "source": [
- "## Adaptive Runge-Kutta Methods\n"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 19,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "name": "stdout",
- "output_type": "stream",
- "text": [
- "'ode23' is a function from the file /home/ryan/octave/odepkg-0.8.5/ode23.m\n",
- "\n",
- " -- Function File: [] = ode23 (@FUN, SLOT, INIT, [OPT], [PAR1, PAR2,\n",
- " ...])\n",
- " -- Command: [SOL] = ode23 (@FUN, SLOT, INIT, [OPT], [PAR1, PAR2, ...])\n",
- " -- Command: [T, Y, [XE, YE, IE]] = ode23 (@FUN, SLOT, INIT, [OPT],\n",
- " [PAR1, PAR2, ...])\n",
- "\n",
- " This function file can be used to solve a set of non-stiff ordinary\n",
- " differential equations (non-stiff ODEs) or non-stiff differential\n",
- " algebraic equations (non-stiff DAEs) with the well known explicit\n",
- " Runge-Kutta method of order (2,3).\n",
- "\n",
- " If this function is called with no return argument then plot the\n",
- " solution over time in a figure window while solving the set of ODEs\n",
- " that are defined in a function and specified by the function handle\n",
- " @FUN. The second input argument SLOT is a double vector that\n",
- " defines the time slot, INIT is a double vector that defines the\n",
- " initial values of the states, OPT can optionally be a structure\n",
- " array that keeps the options created with the command 'odeset' and\n",
- " PAR1, PAR2, ... can optionally be other input arguments of any type\n",
- " that have to be passed to the function defined by @FUN.\n",
- "\n",
- " If this function is called with one return argument then return the\n",
- " solution SOL of type structure array after solving the set of ODEs.\n",
- " The solution SOL has the fields X of type double column vector for\n",
- " the steps chosen by the solver, Y of type double column vector for\n",
- " the solutions at each time step of X, SOLVER of type string for the\n",
- " solver name and optionally the extended time stamp information XE,\n",
- " the extended solution information YE and the extended index\n",
- " information IE all of type double column vector that keep the\n",
- " informations of the event function if an event function handle is\n",
- " set in the option argument OPT.\n",
- "\n",
- " If this function is called with more than one return argument then\n",
- " return the time stamps T, the solution values Y and optionally the\n",
- " extended time stamp information XE, the extended solution\n",
- " information YE and the extended index information IE all of type\n",
- " double column vector.\n",
- "\n",
- " For example, solve an anonymous implementation of the Van der Pol\n",
- " equation\n",
- "\n",
- " fvdb = @(vt,vy) [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];\n",
- "\n",
- " vopt = odeset (\"RelTol\", 1e-3, \"AbsTol\", 1e-3, \\\n",
- " \"NormControl\", \"on\", \"OutputFcn\", @odeplot);\n",
- " ode23 (fvdb, [0 20], [2 0], vopt);\n",
- "\n",
- "See also: odepkg.\n",
- "\n",
- "Additional help for built-in functions and operators is\n",
- "available in the online version of the manual. Use the command\n",
- "'doc ' to search the manual index.\n",
- "\n",
- "Help and information about Octave is also available on the WWW\n",
- "at http://www.octave.org and via the help@octave.org\n",
- "mailing list.\n"
- ]
- }
- ],
- "source": [
- "help ode23"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 20,
- "metadata": {
- "collapsed": false,
- "slideshow": {
- "slide_type": "slide"
- }
- },
- "outputs": [
- {
- "data": {
- "image/svg+xml": [
- "