diff --git a/README.html b/README.html
new file mode 100644
index 0000000..cbdeaa6
--- /dev/null
+++ b/README.html
@@ -0,0 +1,25 @@
+
+
+
+
+
+
+ README
+
+
+
+
+Computational Mechanics 03- Initial Value Problems
+
+
Warning: Under Construction!!
+
+
+
diff --git a/README.md b/README.md
index 1d15de0..b26bad5 100644
--- a/README.md
+++ b/README.md
@@ -1 +1,4 @@
-# Computational Mechanics 03
+# Computational Mechanics 03- Initial Value Problems
+
+
diff --git a/data/.ipynb_checkpoints/therocket-checkpoint.txt b/data/.ipynb_checkpoints/therocket-checkpoint.txt
new file mode 100644
index 0000000..444c7a3
--- /dev/null
+++ b/data/.ipynb_checkpoints/therocket-checkpoint.txt
@@ -0,0 +1,151 @@
+ 0.0000000e+000 2.7316440e-001
+ 1.0000000e-001 1.4411079e+000
+ 2.0000000e-001 2.6693138e+000
+ 3.0000000e-001 4.2383806e+000
+ 4.0000000e-001 5.6499504e+000
+ 5.0000000e-001 6.9489214e+000
+ 6.0000000e-001 8.4234922e+000
+ 7.0000000e-001 9.7402175e+000
+ 8.0000000e-001 1.1202077e+001
+ 9.0000000e-001 1.2641944e+001
+ 1.0000000e+000 1.4119110e+001
+ 1.1000000e+000 1.4143788e+001
+ 1.2000000e+000 1.3868911e+001
+ 1.3000000e+000 1.4228543e+001
+ 1.4000000e+000 1.4349276e+001
+ 1.5000000e+000 1.4332521e+001
+ 1.6000000e+000 1.4485844e+001
+ 1.7000000e+000 1.4544145e+001
+ 1.8000000e+000 1.4663928e+001
+ 1.9000000e+000 1.4723150e+001
+ 2.0000000e+000 1.4762523e+001
+ 2.1000000e+000 1.4567879e+001
+ 2.2000000e+000 1.4700960e+001
+ 2.3000000e+000 1.4935190e+001
+ 2.4000000e+000 1.4835846e+001
+ 2.5000000e+000 1.4939327e+001
+ 2.6000000e+000 1.5135346e+001
+ 2.7000000e+000 1.5135338e+001
+ 2.8000000e+000 1.5306380e+001
+ 2.9000000e+000 1.5132562e+001
+ 3.0000000e+000 1.5381284e+001
+ 3.1000000e+000 1.5236603e+001
+ 3.2000000e+000 1.5322400e+001
+ 3.3000000e+000 1.5562711e+001
+ 3.4000000e+000 1.5585964e+001
+ 3.5000000e+000 1.5553633e+001
+ 3.6000000e+000 1.5704080e+001
+ 3.7000000e+000 1.5741746e+001
+ 3.8000000e+000 1.5777032e+001
+ 3.9000000e+000 1.5958193e+001
+ 4.0000000e+000 1.5851034e+001
+ 4.1000000e+000 1.4431351e+001
+ 4.2000000e+000 1.2597492e+001
+ 4.3000000e+000 1.1252899e+001
+ 4.4000000e+000 9.6343471e+000
+ 4.5000000e+000 8.0758193e+000
+ 4.6000000e+000 6.3308060e+000
+ 4.7000000e+000 4.8680179e+000
+ 4.8000000e+000 3.0927459e+000
+ 4.9000000e+000 1.6899772e+000
+ 5.0000000e+000 -2.1230924e-001
+ 5.1000000e+000 -9.7152877e-001
+ 5.2000000e+000 -2.0733323e+000
+ 5.3000000e+000 -3.0773376e+000
+ 5.4000000e+000 -3.9848158e+000
+ 5.5000000e+000 -5.0336843e+000
+ 5.6000000e+000 -5.9029239e+000
+ 5.7000000e+000 -7.0107236e+000
+ 5.8000000e+000 -7.8986508e+000
+ 5.9000000e+000 -9.0475347e+000
+ 6.0000000e+000 -9.9931052e+000
+ 6.1000000e+000 -9.8490297e+000
+ 6.2000000e+000 -9.6661452e+000
+ 6.3000000e+000 -9.6046216e+000
+ 6.4000000e+000 -9.5843230e+000
+ 6.5000000e+000 -9.5816218e+000
+ 6.6000000e+000 -9.4019202e+000
+ 6.7000000e+000 -9.1890537e+000
+ 6.8000000e+000 -9.2108833e+000
+ 6.9000000e+000 -8.9708582e+000
+ 7.0000000e+000 -8.7781810e+000
+ 7.1000000e+000 -8.7532819e+000
+ 7.2000000e+000 -8.6501691e+000
+ 7.3000000e+000 -8.5149324e+000
+ 7.4000000e+000 -8.3228463e+000
+ 7.5000000e+000 -8.1884909e+000
+ 7.6000000e+000 -8.3247360e+000
+ 7.7000000e+000 -8.0905693e+000
+ 7.8000000e+000 -7.9411118e+000
+ 7.9000000e+000 -7.9152913e+000
+ 8.0000000e+000 -7.5282460e+000
+ 8.1000000e+000 -7.5810719e+000
+ 8.2000000e+000 -7.6406509e+000
+ 8.3000000e+000 -7.3632566e+000
+ 8.4000000e+000 -7.2633092e+000
+ 8.5000000e+000 -7.1462284e+000
+ 8.6000000e+000 -7.2824020e+000
+ 8.7000000e+000 -6.8462979e+000
+ 8.8000000e+000 -7.0498736e+000
+ 8.9000000e+000 -6.6668252e+000
+ 9.0000000e+000 -6.7776371e+000
+ 9.1000000e+000 -6.5170086e+000
+ 9.2000000e+000 -6.3479214e+000
+ 9.3000000e+000 -6.2515036e+000
+ 9.4000000e+000 -6.2185173e+000
+ 9.5000000e+000 -6.2037124e+000
+ 9.6000000e+000 -6.0111919e+000
+ 9.7000000e+000 -5.9691919e+000
+ 9.8000000e+000 -5.9442784e+000
+ 9.9000000e+000 -5.7568067e+000
+ 1.0000000e+001 -5.4967206e+000
+ 1.0100000e+001 -5.3890285e+000
+ 1.0200000e+001 -5.3748506e+000
+ 1.0300000e+001 -5.2160427e+000
+ 1.0400000e+001 -5.0653680e+000
+ 1.0500000e+001 -4.9800986e+000
+ 1.0600000e+001 -4.8631331e+000
+ 1.0700000e+001 -4.5697048e+000
+ 1.0800000e+001 -4.8943903e+000
+ 1.0900000e+001 -4.5216533e+000
+ 1.1000000e+001 -4.4154551e+000
+ 1.1100000e+001 -4.2671072e+000
+ 1.1200000e+001 -4.2803082e+000
+ 1.1300000e+001 -4.0223359e+000
+ 1.1400000e+001 -3.9828129e+000
+ 1.1500000e+001 -3.8040067e+000
+ 1.1600000e+001 -3.6814009e+000
+ 1.1700000e+001 -3.5344749e+000
+ 1.1800000e+001 -3.5619900e+000
+ 1.1900000e+001 -3.3127391e+000
+ 1.2000000e+001 -3.3105316e+000
+ 1.2100000e+001 -3.3651859e+000
+ 1.2200000e+001 -3.1260812e+000
+ 1.2300000e+001 -3.0504968e+000
+ 1.2400000e+001 -3.0618030e+000
+ 1.2500000e+001 -2.8195250e+000
+ 1.2600000e+001 -2.7281635e+000
+ 1.2700000e+001 -2.4834779e+000
+ 1.2800000e+001 -2.4105081e+000
+ 1.2900000e+001 -2.2450488e+000
+ 1.3000000e+001 -2.1937977e+000
+ 1.3100000e+001 -2.1256652e+000
+ 1.3200000e+001 -2.0089646e+000
+ 1.3300000e+001 -1.8599728e+000
+ 1.3400000e+001 -1.6612947e+000
+ 1.3500000e+001 -1.5860937e+000
+ 1.3600000e+001 -1.6911199e+000
+ 1.3700000e+001 -1.4323552e+000
+ 1.3800000e+001 -1.3555511e+000
+ 1.3900000e+001 -1.1650490e+000
+ 1.4000000e+001 -1.1411252e+000
+ 1.4100000e+001 -8.8657232e-001
+ 1.4200000e+001 -9.0682453e-001
+ 1.4300000e+001 -9.2448445e-001
+ 1.4400000e+001 -5.2713209e-001
+ 1.4500000e+001 -5.1147199e-001
+ 1.4600000e+001 -3.8790604e-001
+ 1.4700000e+001 -4.0269559e-001
+ 1.4800000e+001 -1.3883530e-001
+ 1.4900000e+001 -3.3484895e-001
+ 1.5000000e+001 1.0976441e-001
diff --git a/notebooks/1_Catch_Motion.ipynb b/notebooks/.ipynb_checkpoints/01_Catch_Motion-checkpoint.ipynb
similarity index 99%
rename from notebooks/1_Catch_Motion.ipynb
rename to notebooks/.ipynb_checkpoints/01_Catch_Motion-checkpoint.ipynb
index 528aa33..8d8ca10 100644
--- a/notebooks/1_Catch_Motion.ipynb
+++ b/notebooks/.ipynb_checkpoints/01_Catch_Motion-checkpoint.ipynb
@@ -4,7 +4,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2017 L.A. Barba, N.C. Clementi"
+ "###### Modified under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2020 R.C. Cooper"
]
},
{
@@ -13,9 +13,9 @@
"source": [
"# Catch things in motion\n",
"\n",
- "This module of the _Engineering Computations_ course is our launching pad to investigate _change_, _motion_, _dynamics_, using computational thinking, Python, and Jupyter.\n",
+ "This module of the _Computational Mechanics_ course is our launching pad to investigate _change_, _motion_, and _dynamics_, using computational thinking, Python, and Jupyter.\n",
"\n",
- "The foundation of physics and engineering is the subject of **mechanics**: how things move around, when pushed around. Or pulled… in the beginning of the history of mechanics, Galileo and Newton seeked to understand how and why objects fall under the pull of gravity.\n",
+ "The foundation of physics and engineering is the subject of **mechanics**: how things move around, when pushed around. Or pulled... in the beginning of the history of mechanics, Galileo and Newton sought to understand how and why objects fall under the pull of gravity.\n",
"\n",
"This first lesson will explore motion by analyzing images and video, to learn about velocity and acceleration."
]
@@ -33,7 +33,7 @@
},
{
"cell_type": "code",
- "execution_count": 1,
+ "execution_count": 13,
"metadata": {},
"outputs": [
{
@@ -51,7 +51,7 @@
" "
],
"text/plain": [
- ""
+ ""
]
},
"metadata": {},
@@ -111,7 +111,7 @@
},
{
"cell_type": "code",
- "execution_count": 4,
+ "execution_count": 28,
"metadata": {},
"outputs": [
{
@@ -120,7 +120,7 @@
"imageio.plugins.ffmpeg.FfmpegFormat.Reader"
]
},
- "execution_count": 4,
+ "execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
@@ -131,27 +131,7 @@
},
{
"cell_type": "code",
- "execution_count": 5,
- "metadata": {},
- "outputs": [
- {
- "data": {
- "text/plain": [
- "1233"
- ]
- },
- "execution_count": 5,
- "metadata": {},
- "output_type": "execute_result"
- }
- ],
- "source": [
- "len(reader)"
- ]
- },
- {
- "cell_type": "code",
- "execution_count": 6,
+ "execution_count": 31,
"metadata": {},
"outputs": [
{
@@ -173,7 +153,7 @@
"source": [
"##### Note:\n",
"\n",
- "You may get this error after calling `get_reader()`:\n",
+ "You may get this error after calling `get_reader()` if you're not running in the class Jupyter server:\n",
" \n",
"```\n",
"NeedDownloadError: Need ffmpeg exe. You can obtain it with either:\n",
diff --git a/notebooks/.ipynb_checkpoints/2_Step_Future-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/02_Step_Future-checkpoint.ipynb
similarity index 100%
rename from notebooks/.ipynb_checkpoints/2_Step_Future-checkpoint.ipynb
rename to notebooks/.ipynb_checkpoints/02_Step_Future-checkpoint.ipynb
diff --git a/notebooks/3_Get_Oscillations.ipynb b/notebooks/.ipynb_checkpoints/03_Get_Oscillations-checkpoint.ipynb
similarity index 100%
rename from notebooks/3_Get_Oscillations.ipynb
rename to notebooks/.ipynb_checkpoints/03_Get_Oscillations-checkpoint.ipynb
diff --git a/notebooks/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb
new file mode 100644
index 0000000..aa2b764
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/06_roots-1-checkpoint.ipynb
@@ -0,0 +1,892 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "metadata": {
+ "collapsed": true,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Roots of Nonlinear functions\n",
+ "## Bracketing ch. 5"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "## Not always possible to solve for a given variable. \n",
+ "\n",
+ "### Freefall example:\n",
+ "If an object, with drag coefficient of 0.25 kg/m reaches a velocity of 36 m/s after 4 seconds of freefalling, what is its mass?\n",
+ "\n",
+ "$v(t)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)$\n",
+ "\n",
+ "Cannot solve for m \n",
+ "\n",
+ "Instead, solve the problem by creating a new function f(m) where\n",
+ "\n",
+ "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$. \n",
+ "\n",
+ "When f(m) = 0, we have solved for m in terms of the other variables (e.g. for a given time, velocity, drag coefficient and acceleration due to gravity)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "setdefaults\n",
+ "g=9.81; % acceleration due to gravity\n",
+ "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n",
+ "c_d=0.25; % drag coefficient\n",
+ "t=4; % at time = 4 seconds\n",
+ "v=36; % speed must be 36 m/s\n",
+ "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n",
+ "\n",
+ "plot(m,f_m(m),m,zeros(length(m),1))\n",
+ "axis([45 200 -5 1])"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 0.045626\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "f_m(145)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "Brute force method is plot f_m vs m and with smaller and smaller steps until f_m ~ 0\n",
+ "\n",
+ "Better methods are the \n",
+ "1. Bracketing methods\n",
+ "2. Open methods\n",
+ "\n",
+ "Both need an initial guess. \n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Incremental method (Brute force)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "You know that for one value, m_lower, f_m is negative and for another value, m_upper, f_m is positive. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "^C\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "function xb = incsearch(func,xmin,xmax,ns)\n",
+ "% incsearch: incremental search root locator\n",
+ "% xb = incsearch(func,xmin,xmax,ns):\n",
+ "% finds brackets of x that contain sign changes\n",
+ "% of a function on an interval\n",
+ "% input:\n",
+ "% func = name of function\n",
+ "% xmin, xmax = endpoints of interval\n",
+ "% ns = number of subintervals (default = 50)\n",
+ "% output:\n",
+ "% xb(k,1) is the lower bound of the kth sign change\n",
+ "% xb(k,2) is the upper bound of the kth sign change\n",
+ "% If no brackets found, xb = [].\n",
+ "if nargin < 3, error('at least 3 arguments required'), end\n",
+ "if nargin < 4, ns = 50; end %if ns blank set to 50\n",
+ "% Incremental search\n",
+ "x = linspace(xmin,xmax,ns);\n",
+ "f = func(x);\n",
+ "nb = 0; xb = []; %xb is null unless sign change detected\n",
+ "for k = 1:length(x)-1\n",
+ " if sign(f(k)) ~= sign(f(k+1)) %check for sign change\n",
+ " nb = nb + 1;\n",
+ " xb(nb,1) = x(k);\n",
+ " xb(nb,2) = x(k+1);\n",
+ " end\n",
+ "end\n",
+ "if isempty(xb) %display that no brackets were found\n",
+ " fprintf('no brackets found\\n')\n",
+ " fprintf('check interval or increase ns\\n')\n",
+ "else\n",
+ " fprintf('number of brackets: %i\\n',nb) %display number of brackets\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "number of brackets: 1\n",
+ "ans =\n",
+ "\n",
+ " 141.84 144.90\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "incsearch(f_m,50,200)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "number of brackets: 1\n",
+ "ans =\n",
+ "\n",
+ " 142.73 142.83\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "incsearch(f_m,140, 150,100)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Bisection method"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Divide interval in half until error is reduced to some level\n",
+ "\n",
+ "in previous example of freefall, choose x_l=50, x_u=200\n",
+ "\n",
+ "x_r = (50+200)/2 = 125\n",
+ "\n",
+ "f_m(125) = -0.408\n",
+ "\n",
+ "x_r= (125+200)/2 = 162.5\n",
+ "\n",
+ "f_m(162.5) = 0.3594\n",
+ "\n",
+ "x_r = (125+162.5)/2=143.75\n",
+ "\n",
+ "f_m(143.75)= 0.0206"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 134.38\n",
+ "interval left f(x_l)= -0.4,f(x_r)= -0.2\n",
+ "interval right f(x_r)= -0.2,f(x_u)= 0.0\n",
+ "ans = -0.18060\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_l=125; x_u=143.75;\n",
+ "x_r=(x_l+x_u)/2\n",
+ "fprintf('interval left f(x_l)= %1.1f,f(x_r)= %1.1f\\n',f_m(x_l),f_m(x_r))\n",
+ "fprintf('interval right f(x_r)= %1.1f,f(x_u)= %1.1f\\n',f_m(x_r),f_m(x_u))\n",
+ "f_m(x_r)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Bisect Function"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Much better root locator, with 4 iterations, our function is already close to zero\n",
+ "\n",
+ "Automate this with a function:\n",
+ "`bisect.m`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "^C\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)\n",
+ "% bisect: root location zeroes\n",
+ "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
+ "% uses bisection method to find the root of func\n",
+ "% input:\n",
+ "% func = name of function\n",
+ "% xl, xu = lower and upper guesses\n",
+ "% es = desired relative error (default = 0.0001%)\n",
+ "% maxit = maximum allowable iterations (default = 50)\n",
+ "% p1,p2,... = additional parameters used by func\n",
+ "% output:\n",
+ "% root = real root\n",
+ "% fx = function value at root\n",
+ "% ea = approximate relative error (%)\n",
+ "% iter = number of iterations\n",
+ "if nargin<3,error('at least 3 input arguments required'),end\n",
+ "test = func(xl,varargin{:})*func(xu,varargin{:});\n",
+ "if test>0,error('no sign change'),end\n",
+ "if nargin<4|isempty(es), es=0.0001;end\n",
+ "if nargin<5|isempty(maxit), maxit=50;end\n",
+ "iter = 0; xr = xl; ea = 100;\n",
+ "while (1)\n",
+ " xrold = xr;\n",
+ " xr = (xl + xu)/2;\n",
+ " iter = iter + 1;\n",
+ " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n",
+ " test = func(xl,varargin{:})*func(xr,varargin{:});\n",
+ " if test < 0\n",
+ " xu = xr;\n",
+ " elseif test > 0\n",
+ " xl = xr;\n",
+ " else\n",
+ " ea = 0;\n",
+ " end\n",
+ " if ea <= es | iter >= maxit,break,end\n",
+ "end\n",
+ "root = xr; fx = func(xr, varargin{:});"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Mass_at_36ms = 142.74\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "Mass_at_36ms=bisect(f_m,50,200)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Thanks"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## False position (linear interpolation)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Rather than bisecting each bracket (1/2 each time) we can calculate the slope between the two points and update the xr position in this manner\n",
+ "\n",
+ "$ x_{r} = x_{u} - \\frac{f(x_{u})(x_{l}-x_{u})}{f(x_{l})-f(x_{u})}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {
+ "collapsed": false,
+ "scrolled": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "%xl=50; xu=200; \n",
+ "xu=xr;\n",
+ "xl=xl;\n",
+ "xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu));\n",
+ "plot(m,f_m(m),xl,f_m(xl),'s',xu,f_m(xu),'s',xr,0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## False Position"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Much better root locator, with 4 iterations, our function is already close to zero\n",
+ "\n",
+ "Automate this with a function:\n",
+ "`falsepos.m`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "^C\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "function [root,fx,ea,iter]=falsepos(func,xl,xu,es,maxit,varargin)\n",
+ "% falsepos: root location zeroes\n",
+ "% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):\n",
+ "% uses false position method to find the root of func\n",
+ "% input:\n",
+ "% func = name of function\n",
+ "% xl, xu = lower and upper guesses\n",
+ "% es = desired relative error (default = 0.0001%)\n",
+ "% maxit = maximum allowable iterations (default = 50)\n",
+ "% p1,p2,... = additional parameters used by func\n",
+ "% output:\n",
+ "% root = real root\n",
+ "% fx = function value at root\n",
+ "% ea = approximate relative error (%)\n",
+ "% iter = number of iterations\n",
+ "if nargin<3,error('at least 3 input arguments required'),end\n",
+ "test = func(xl,varargin{:})*func(xu,varargin{:});\n",
+ "if test>0,error('no sign change'),end\n",
+ "if nargin<4|isempty(es), es=0.0001;end\n",
+ "if nargin<5|isempty(maxit), maxit=50;end\n",
+ "iter = 0; xr = xl; ea = 100;\n",
+ "while (1)\n",
+ " xrold = xr;\n",
+ " % xr = (xl + xu)/2; % bisect method\n",
+ " xr=xu - (f_m(xu)*(xl-xu))/(f_m(xl)-f_m(xu)); % false position method\n",
+ " iter = iter + 1;\n",
+ " if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end\n",
+ " test = func(xl,varargin{:})*func(xr,varargin{:});\n",
+ " if test < 0\n",
+ " xu = xr;\n",
+ " elseif test > 0\n",
+ " xl = xr;\n",
+ " else\n",
+ " ea = 0;\n",
+ " end\n",
+ " if ea <= es | iter >= maxit,break,end\n",
+ "end\n",
+ "root = xr; fx = func(xr, varargin{:});"
+ ]
+ }
+ ],
+ "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/notebooks/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb
new file mode 100644
index 0000000..4a2ba51
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/07_roots_and_optimization-2-checkpoint.ipynb
@@ -0,0 +1,1525 @@
+{
+ "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": [
+ "# Roots: Open methods\n",
+ "## Newton-Raphson"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "First-order approximation for the location of the root (i.e. assume the slope at the given point is constant, what is the solution when f(x)=0)\n",
+ "\n",
+ "$f'(x_{i})=\\frac{f(x_{i})-0}{x_{i}-x_{i+1}}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "Use Newton-Raphson to find solution when $e^{-x}=x$"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "source": [
+ "Find x when $e^{-x}=x$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.50000\n",
+ "error_approx = 1\n"
+ ]
+ }
+ ],
+ "source": [
+ "f= @(x) exp(-x)-x;\n",
+ "df= @(x) -exp(-x)-1;\n",
+ "\n",
+ "x_i= 0;\n",
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56714\n",
+ "error_approx = 0.0014673\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56714\n",
+ "error_approx = 2.2106e-07\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "x_r = 0.56714\n",
+ "error_approx = 5.0897e-15\n"
+ ]
+ }
+ ],
+ "source": [
+ "x_r = x_i-f(x_i)/df(x_i)\n",
+ "error_approx = abs((x_r-x_i)/x_r)\n",
+ "x_i=x_r;"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Create function `newtraph.m`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "^C\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "function [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,varargin)\n",
+ "% newtraph: Newton-Raphson root location zeroes\n",
+ "% [root,ea,iter]=newtraph(func,dfunc,xr,es,maxit,p1,p2,...):\n",
+ "% uses Newton-Raphson method to find the root of func\n",
+ "% input:\n",
+ "% func = name of function\n",
+ "% dfunc = name of derivative of function\n",
+ "% xr = initial guess\n",
+ "% es = desired relative error (default = 0.0001%)\n",
+ "% maxit = maximum allowable iterations (default = 50)\n",
+ "% p1,p2,... = additional parameters used by function\n",
+ "% output:\n",
+ "% root = real root\n",
+ "% ea = approximate relative error (%)\n",
+ "% iter = number of iterations\n",
+ "if nargin<3,error('at least 3 input arguments required'),end\n",
+ "if nargin<4 || isempty(es),es=0.0001;end\n",
+ "if nargin<5 || isempty(maxit),maxit=50;end\n",
+ "iter = 0;\n",
+ "while (1)\n",
+ " xrold = xr;\n",
+ " xr = xr - func(xr)/dfunc(xr);\n",
+ " iter = iter + 1;\n",
+ " if xr ~= 0 \n",
+ " ea = abs((xr - xrold)/xr) * 100;\n",
+ " end\n",
+ " if ea <= es || iter >= maxit, break, end\n",
+ "end\n",
+ "root = xr;"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "In the freefall example, we created a function f(m) that when f(m)=0, then the mass had been chosen such that at t=4 s, the velocity is 36 m/s. \n",
+ "\n",
+ "$f(m)=\\sqrt{\\frac{gm}{c_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-v(t)$.\n",
+ "\n",
+ "to use the Newton-Raphson method, we need the derivative $\\frac{df}{dm}$\n",
+ "\n",
+ "$\\frac{df}{dm}=\\frac{1}{2}\\sqrt{\\frac{g}{mc_{d}}}\\tanh(\\sqrt{\\frac{gc_{d}}{m}}t)-\n",
+ "\\frac{g}{2m}\\mathrm{sech}^{2}(\\sqrt{\\frac{gc_{d}}{m}}t)$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults\n",
+ "g=9.81; % acceleration due to gravity\n",
+ "m=linspace(50, 200,100); % possible values for mass 50 to 200 kg\n",
+ "c_d=0.25; % drag coefficient\n",
+ "t=4; % at time = 4 seconds\n",
+ "v=36; % speed must be 36 m/s\n",
+ "f_m = @(m) sqrt(g*m/c_d).*tanh(sqrt(g*c_d./m)*t)-v; % anonymous function f_m\n",
+ "df_m = @(m) 1/2*sqrt(g./m/c_d).*tanh(sqrt(g*c_d./m)*t)-g/2./m*sech(sqrt(g*c_d./m)*t).^2;"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "root = 142.74\n",
+ "ea = 1.8806e-04\n",
+ "iter = 50\n"
+ ]
+ }
+ ],
+ "source": [
+ "[root,ea,iter]=newtraph(f_m,df_m,50,0.0001)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Thanks"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Secant Methods"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Not always able to evaluate the derivative. Approximation of derivative:\n",
+ "\n",
+ "$f'(x_{i})=\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{f'(x_{i})}$\n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})}{\\frac{f(x_{i-1})-f(x_{i})}{x_{i-1}-x_{i}}}=\n",
+ " x_{i}-\\frac{f(x_{i})(x_{i-1}-x_{i})}{f(x_{i-1})-f(x_{i})}$\n",
+ " \n",
+ "What values should $x_{i}$ and $x_{i-1}$ take?\n",
+ "\n",
+ "To reduce arbitrary selection of variables, use the"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Modified Secant method"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "fragment"
+ }
+ },
+ "source": [
+ "Change the x evaluations to a perturbation $\\delta$. \n",
+ "\n",
+ "$x_{i+1}=x_{i}-\\frac{f(x_{i})(\\delta x_{i})}{f(x_{i}+\\delta x_{i})-f(x_{i})}$"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "root = 142.74\n",
+ "ea = 3.0615e-07\n",
+ "iter = 7\n"
+ ]
+ }
+ ],
+ "source": [
+ "[root,ea,iter]=mod_secant(f_m,1,50,0.00001)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.1185e+04\r\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "car_payments(400,30000,0.05,5,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 17,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Amt_numerical = 5467.0\n",
+ "ans = 3.9755e-04\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "Amt_numerical=mod_secant(@(A) car_payments(A,700000,0.0875,30,0),1e-6,50,0.001)\n",
+ "car_payments(Amt_numerical,700000,0.0875,30,1)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1.9681e+06\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "Amt_numerical*12*30"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Amortization\n",
+ "\n",
+ "Amortization calculation makes the same calculation for the monthly payment amount, A, paying off the principle amount, P, over n pay periods with monthly interest rate, r. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "Amt = 566.14\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "% Amortization calculation\n",
+ "A = @(P,r,n) P*(r*(1+r)^n)./((1+r)^n-1);\n",
+ "Amt=A(30000,0.05/12,5*12)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Matlab's function\n",
+ "\n",
+ "Matlab and Octave combine bracketing and open methods in the `fzero` function. "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "'fzero' is a function from the file /usr/share/octave/4.0.0/m/optimization/fzero.m\n",
+ "\n",
+ " -- Function File: fzero (FUN, X0)\n",
+ " -- Function File: fzero (FUN, X0, OPTIONS)\n",
+ " -- Function File: [X, FVAL, INFO, OUTPUT] = fzero (...)\n",
+ " Find a zero of a univariate function.\n",
+ "\n",
+ " FUN is a function handle, inline function, or string containing the\n",
+ " name of the function to evaluate.\n",
+ "\n",
+ " X0 should be a two-element vector specifying two points which\n",
+ " bracket a zero. In other words, there must be a change in sign of\n",
+ " the function between X0(1) and X0(2). More mathematically, the\n",
+ " following must hold\n",
+ "\n",
+ " sign (FUN(X0(1))) * sign (FUN(X0(2))) <= 0\n",
+ "\n",
+ " If X0 is a single scalar then several nearby and distant values are\n",
+ " probed in an attempt to obtain a valid bracketing. If this is not\n",
+ " successful, the function fails.\n",
+ "\n",
+ " OPTIONS is a structure specifying additional options. Currently,\n",
+ " 'fzero' recognizes these options: \"FunValCheck\", \"OutputFcn\",\n",
+ " \"TolX\", \"MaxIter\", \"MaxFunEvals\". For a description of these\n",
+ " options, see *note optimset: XREFoptimset.\n",
+ "\n",
+ " On exit, the function returns X, the approximate zero point and\n",
+ " FVAL, the function value thereof.\n",
+ "\n",
+ " INFO is an exit flag that can have these values:\n",
+ "\n",
+ " * 1 The algorithm converged to a solution.\n",
+ "\n",
+ " * 0 Maximum number of iterations or function evaluations has\n",
+ " been reached.\n",
+ "\n",
+ " * -1 The algorithm has been terminated from user output\n",
+ " function.\n",
+ "\n",
+ " * -5 The algorithm may have converged to a singular point.\n",
+ "\n",
+ " OUTPUT is a structure containing runtime information about the\n",
+ " 'fzero' algorithm. Fields in the structure are:\n",
+ "\n",
+ " * iterations Number of iterations through loop.\n",
+ "\n",
+ " * nfev Number of function evaluations.\n",
+ "\n",
+ " * bracketx A two-element vector with the final bracketing of the\n",
+ " zero along the x-axis.\n",
+ "\n",
+ " * brackety A two-element vector with the final bracketing of the\n",
+ " zero along the y-axis.\n",
+ "\n",
+ " See also: optimset, fsolve.\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 fzero"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 563.79\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "fzero(@(A) car_payments(A,30000,0.05,5,0),500)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "## Comparison of Solvers\n",
+ "\n",
+ "It's helpful to compare to the convergence of different routines to see how quickly you find a solution. \n",
+ "\n",
+ "Comparing the freefall example\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 23,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: called from\n",
+ " __line__ at line 120 column 16\n",
+ " line at line 56 column 8\n",
+ " __plt__>__plt2vv__ at line 500 column 10\n",
+ " __plt__>__plt2__ at line 246 column 14\n",
+ " __plt__ at line 133 column 15\n",
+ " semilogy at line 60 column 10\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=20;\n",
+ "iterations = linspace(1,400,N);\n",
+ "ea_nr=zeros(1,N); % appr error Newton-Raphson\n",
+ "ea_ms=zeros(1,N); % appr error Modified Secant\n",
+ "ea_fp=zeros(1,N); % appr error false point method\n",
+ "ea_bs=zeros(1,N); % appr error bisect method\n",
+ "for i=1:length(iterations)\n",
+ " [root_nr,ea_nr(i),iter_nr]=newtraph(f_m,df_m,300,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f_m,1e-6,300,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f_m,1,300,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f_m,1,300,0,iterations(i));\n",
+ "end\n",
+ "\n",
+ "setdefaults\n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 22,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "subslide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_nr =\n",
+ "\n",
+ " Columns 1 through 8:\n",
+ "\n",
+ " 6.36591 0.06436 0.00052 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 9 through 16:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ " Columns 17 through 20:\n",
+ "\n",
+ " 0.00000 0.00000 0.00000 0.00000\n",
+ "\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_nr"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 26,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: called from\n",
+ " __line__ at line 120 column 16\n",
+ " line at line 56 column 8\n",
+ " __plt__>__plt2vv__ at line 500 column 10\n",
+ " __plt__>__plt2__ at line 246 column 14\n",
+ " __plt__ at line 133 column 15\n",
+ " semilogy at line 60 column 10\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n",
+ "warning: axis: omitting non-positive data in log plot\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "N=20;\n",
+ "f= @(x) x^10-1;\n",
+ "df=@(x) 10*x^9;\n",
+ "iterations = linspace(1,50,N);\n",
+ "ea_nr=zeros(1,N); % appr error Newton-Raphson\n",
+ "ea_ms=zeros(1,N); % appr error Modified Secant\n",
+ "ea_fp=zeros(1,N); % appr error false point method\n",
+ "ea_bs=zeros(1,N); % appr error bisect method\n",
+ "for i=1:length(iterations)\n",
+ " [root_nr,ea_nr(i),iter_nr]=newtraph(f,df,0.5,0,iterations(i));\n",
+ " [root_ms,ea_ms(i),iter_ms]=mod_secant(f,1e-6,0.5,0,iterations(i));\n",
+ " [root_fp,ea_fp(i),iter_fp]=falsepos(f,0,5,0,iterations(i));\n",
+ " [root_bs,ea_bs(i),iter_bs]=bisect(f,0,5,0,iterations(i));\n",
+ "end\n",
+ " \n",
+ "semilogy(iterations,abs(ea_nr),iterations,abs(ea_ms),iterations,abs(ea_fp),iterations,abs(ea_bs))\n",
+ "legend('newton-raphson','mod-secant','false point','bisection')"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 27,
+ "metadata": {
+ "collapsed": false,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ea_nr =\n",
+ "\n",
+ " Columns 1 through 7:\n",
+ "\n",
+ " 99.03195 11.11111 11.11111 11.11111 11.11111 11.11111 11.11111\n",
+ "\n",
+ " Columns 8 through 14:\n",
+ "\n",
+ " 11.11111 11.11111 11.11111 11.11109 11.11052 11.10624 10.99684\n",
+ "\n",
+ " Columns 15 through 20:\n",
+ "\n",
+ " 8.76956 2.12993 0.00000 0.00000 0.00000 0.00000\n",
+ "\n",
+ "ans = 16.208\n"
+ ]
+ }
+ ],
+ "source": [
+ "ea_nr\n",
+ "newtraph(f,df,0.5,0,12)"
+ ]
+ }
+ ],
+ "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/notebooks/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb
new file mode 100644
index 0000000..e67cb49
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/18_initial_value_ode-checkpoint.ipynb
@@ -0,0 +1,794 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "skip"
+ }
+ },
+ "outputs": [],
+ "source": [
+ "setdefaults"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Initial Value Problems (ODEs)\n",
+ "*Ch. 22*"
+ ]
+ },
+ {
+ "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": 11,
+ "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",
+ " xlabel('time (s)')\n",
+ " ylabel('velocity (m/s)')\n",
+ "end"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "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",
+ " 6.7 | 23.01\n",
+ " 6.7 | 23.01\n",
+ " 13.3 | 36.09\n",
+ " 13.3 | 36.09\n",
+ " 20.0 | 24.90\n",
+ " 20.0 | 24.90\n",
+ "\n",
+ "O(h^2)=44.44\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[v_an,v_t,t]=freefall(4,20);\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": 19,
+ "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": 20,
+ "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": 21,
+ "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": 22,
+ "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": 24,
+ "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": "markdown",
+ "metadata": {
+ "collapsed": true,
+ "slideshow": {
+ "slide_type": "slide"
+ }
+ },
+ "source": [
+ "# Thanks"
+ ]
+ }
+ ],
+ "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/notebooks/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb b/notebooks/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb
new file mode 100644
index 0000000..361a7bc
--- /dev/null
+++ b/notebooks/.ipynb_checkpoints/19_nonlinear_ivp-checkpoint.ipynb
@@ -0,0 +1,3482 @@
+{
+ "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": [
+ "