diff --git a/lecture_02/.ipynb_checkpoints/lecture_02-checkpoint.ipynb b/lecture_02/.ipynb_checkpoints/lecture_02-checkpoint.ipynb
new file mode 100644
index 0000000..e34d82f
--- /dev/null
+++ b/lecture_02/.ipynb_checkpoints/lecture_02-checkpoint.ipynb
@@ -0,0 +1,594 @@
+{
+ "cells": [
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Solution to Form #1"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 21,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 14\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "[1,2,3]*[1;2;3]"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "[1,2,3]*[1;2;3]=?"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# The first source of error is roundoff error\n",
+ "## Just storing a number in a computer requires rounding"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "realmax = 1.79769313486231570815e+308\n",
+ "realmin = 2.22507385850720138309e-308\n",
+ "maximum relative error = 2.22044604925031308085e-16\n"
+ ]
+ }
+ ],
+ "source": [
+ "fprintf('realmax = %1.20e\\n',realmax)\n",
+ "fprintf('realmin = %1.20e\\n',realmin)\n",
+ "fprintf('maximum relative error = %1.20e\\n',eps)\n"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "ans = 1\r\n"
+ ]
+ }
+ ],
+ "source": [
+ "s=1;\n",
+ "for i=1:1000\n",
+ " s=s+eps/10;\n",
+ "end\n",
+ "s==1"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Freefall Model (revisited)\n",
+ "## Octave solution (will run same on Matlab)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Set default values in Octave for linewidth and text size"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "%plot --format svg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 15,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "set (0, \"defaultaxesfontname\", \"Helvetica\")\n",
+ "set (0, \"defaultaxesfontsize\", 18)\n",
+ "set (0, \"defaulttextfontname\", \"Helvetica\")\n",
+ "set (0, \"defaulttextfontsize\", 18) \n",
+ "set (0, \"defaultlinelinewidth\", 4)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Define time from 0 to 12 seconds with `N` timesteps \n",
+ "function defined as `freefall`"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [],
+ "source": [
+ "function [v_analytical,v_terminal,t]=freefall(N)\n",
+ " t=linspace(0,12,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)|vel analytical (m/s)|vel numerical (m/s)\\n')\n",
+ " fprintf('-----------------------------------------------\\n')\n",
+ " M=[t(indices),v_analytical(indices),v_numerical(indices)];\n",
+ " fprintf('%7.1f | %18.2f | %15.2f\\n',M(:,1:3)');\n",
+ " plot(t,v_analytical,'-',t,v_numerical,'o-')\n",
+ "end\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 19,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "time (s)|vel analytical (m/s)|vel numerical (m/s)\n",
+ "-----------------------------------------------\n",
+ " 0.0 | 0.00 | 0.00\n",
+ " 2.0 | 18.76 | 18.82\n",
+ " 4.0 | 32.64 | 32.80\n",
+ " 6.1 | 40.79 | 40.97\n",
+ " 8.0 | 44.80 | 44.94\n",
+ " 10.0 | 46.84 | 46.93\n",
+ " 12.0 | 47.77 | 47.82\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[v_analytical,v_terminal,t]=freefall(120);"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "# Types of error\n",
+ "## Freefall is example of \"truncation error\"\n",
+ "### Truncation error results from approximating exact mathematical procedure"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "We approximated the derivative as $\\delta v/\\delta t\\approx\\Delta v/\\Delta t$\n",
+ "\n",
+ "Can reduce error by decreasing step size -> $\\Delta t$=`delta_time`"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "## Another example of truncation error is a Taylor series (or Maclaurin if centered at a=0)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Taylor series:\n",
+ "$f(x)=f(a)+f'(a)(x-a)+\\frac{f''(a)}{2!}(x-a)^{2}+\\frac{f'''(a)}{3!}(x-a)^{3}+...$\n",
+ "\n",
+ "We can approximate the next value in a function by adding Taylor series terms:\n",
+ "\n",
+ "|Approximation | formula |\n",
+ "|---|-------------------------|\n",
+ "|$0^{th}$-order | $f(x_{i+1})=f(x_{i})+R_{1}$ |\n",
+ "|$1^{st}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+R_{2}$ |\n",
+ "|$2^{nd}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+\\frac{f''(x_{i})}{2!}h^{2}+R_{3}$|\n",
+ "|$n^{th}$-order | $f(x_{i+1})=f(x_{i})+f'(x_{i})h+\\frac{f''(x_{i})}{2!}h^{2}+...\\frac{f^{(n)}}{n!}h^{n}+R_{n}$|\n",
+ "\n",
+ "Where $R_{n}=\\frac{f^{(n+1)}(\\xi)}{(n+1)!}h^{n+1}$ is the error associated with truncating the approximation at order $n$.\n",
+ "\n",
+ "The $n^{th}$-order approximation estimates that the unknown function, $f(x)$, is equal to an $n^{th}$-order polynomial. \n",
+ "\n",
+ "In the Freefall example, we estimated the function with a $1^{st}$-order approximation, so \n",
+ "\n",
+ "$v(t_{i+1})=v(t_{i})+v'(t_{i})(t_{i+1}-t_{i})+R_{1}$\n",
+ "\n",
+ "$v'(t_{i})=\\frac{v(t_{i+1})-v(t_{i})}{t_{i+1}-t_{i}}-\\frac{R_{1}}{t_{i+1}-t_{i}}$\n",
+ "\n",
+ "$\\frac{R_{1}}{t_{i+1}-t_{i}}=\\frac{v''(\\xi)}{2!}(t_{i+1}-t_{i})$\n",
+ "\n",
+ "or the truncation error for a first-order Taylor series approximation is\n",
+ "\n",
+ "$\\frac{R_{1}}{t_{i+1}-t_{i}}=O(\\Delta t)$\n"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "1. digital representation of a number is rarely exact\n",
+ "\n",
+ "2. arithmetic (+,-,/,\\*) causes roundoff error"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 20,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "3.14159265358979311600\n",
+ "3.14159274101257324219\n"
+ ]
+ }
+ ],
+ "source": [
+ "fprintf('%1.20f\\n',double(pi))\n",
+ "fprintf('%1.20f\\n',single(pi))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "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/lecture_21/lecture_21.ipynb b/lecture_21/lecture_21.ipynb
index c2547e6..6d1701d 100644
--- a/lecture_21/lecture_21.ipynb
+++ b/lecture_21/lecture_21.ipynb
@@ -55,7 +55,7 @@
},
{
"cell_type": "code",
- "execution_count": 6,
+ "execution_count": 3,
"metadata": {
"collapsed": false
},
@@ -208,111 +208,111 @@
"\n",
"\n",
"\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\tgnuplot_plot_2a\n",
"\n",
"\t \n",
- "\t\n",
- "\t\n",
+ "\t\n",
"\t\n",
"\n",
@@ -393,7 +393,750 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "# Initial Value Problem"
+ "# 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": 4,
+ "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": 5,
+ "metadata": {
+ "collapsed": false
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "time (s)| error (m/s)\n",
+ "-------------------------\n",
+ " 0.0 | 0.00\n",
+ " 1.7 | 0.64\n",
+ " 3.4 | 2.50\n",
+ " 6.9 | 3.12\n",
+ " 8.6 | 2.10\n",
+ " 10.3 | 1.23\n",
+ " 12.0 | 0.67\n",
+ "\n",
+ "O(h^2)=2.94\n"
+ ]
+ },
+ {
+ "data": {
+ "image/svg+xml": [
+ ""
+ ],
+ "text/plain": [
+ ""
+ ]
+ },
+ "metadata": {},
+ "output_type": "display_data"
+ }
+ ],
+ "source": [
+ "[v_an,v_t,t]=freefall(8,12);\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')"
]
},
{
diff --git a/lecture_21/octave-workspace b/lecture_21/octave-workspace
new file mode 100644
index 0000000..398435f
Binary files /dev/null and b/lecture_21/octave-workspace differ