added HW6 mini project competition

19_nonlinear_ivp/19_nonlinear_ivp.ipynb

@@ -2,7 +2,7 @@
  "cells": [
   {
    "cell_type": "code",
-   "execution_count": 1,
+   "execution_count": null,
    "metadata": {
     "collapsed": true,
     "slideshow": {
@@ -16,7 +16,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": null,
    "metadata": {
     "collapsed": true,
     "slideshow": {
@@ -30,7 +30,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": null,
    "metadata": {
     "collapsed": false,
     "slideshow": {
@@ -51,8 +51,17 @@
    },
    "source": [
     "# Initial Value Problems (continued)\n",
-    "*ch. 23 - Adaptive methods*\n",
-    "\n",
+    "*ch. 23 - Adaptive methods*"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "slideshow": {
+     "slide_type": "slide"
+    }
+   },
+   "source": [
     "## Predator-Prey Models (Chaos Theory)"
    ]
   },
@@ -77,7 +86,11 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "slideshow": {
+     "slide_type": "slide"
+    }
+   },
    "source": [
     "```matlab\n",
     "function yp=predprey(t,y,a,b,c,d)\n",
@@ -93,28 +106,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": null,
    "metadata": {
     "collapsed": false,
     "slideshow": {
      "slide_type": "fragment"
     }
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "ans =\n",
-      "\n",
-      "   12.0000\n",
-      "    5.2000\n",
-      "\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "predprey(1,[20 1],a,b,c,d)"
+    "predprey(1,[20 1],1,1,1,1)"
    ]
   },
   {
@@ -515,6 +516,73 @@
     "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)<tf\n",
+    "    t(n+1) = tf;\n",
+    "    n = n+1;\n",
+    "  end\n",
+    "else\n",
+    "  t = tspan;\n",
+    "end\n",
+    "tt = ti; y(1,:) = y0;\n",
+    "np = 1; tp(np) = tt; yp(np,:) = y(1,:);\n",
+    "i=1;\n",
+    "while(1)\n",
+    "  tend = t(np+1);\n",
+    "  hh = t(np+1) - t(np);\n",
+    "  if hh>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": 43,
@@ -992,13 +1060,13 @@
    ]
   },
   {
-   "cell_type": "code",
-   "execution_count": null,
+   "cell_type": "markdown",
    "metadata": {
     "collapsed": true
    },
-   "outputs": [],
-   "source": []
+   "source": [
+    "# Thanks"
+   ]
   }
  ],
  "metadata": {

HW6/README.md

@@ -41,6 +41,85 @@ b. with Euler's method
 
 c. with Heun's predictor-corrector approach
 
-**5\.** *Coming Monday afternoon*
+### Paper airplane  Mini-design challenge 
+[phugoid
+models in Python](http://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/01_03_PhugoidFullModel.ipynb)
+
+![Image](http://nbviewer.jupyter.org/github/numerical-mooc/numerical-mooc/blob/master/lessons/01_phugoid/figures/glider_forces-lesson3.png)
+
+**5\.** 
+
+Using a phugoid model, write a function and script to analyze the flight of a paper airplane.
+
+
+![eq1](./equations/eq1.png)
+
+Use primes to
+denote the time derivatives and divide through by the weight:
+
+![eq2](./equations/eq2.png)
+
+
+Ratio of lift to weight is known from the trim
+ conditions ![eq3](./equations/eq3.png) and also from the definitions of lift and drag, 
+
+ ![eq4](./equations/eq4.png) 
+
+we see that L/D=C_L/C_D. The system of equations can be re-written:
+
+![eq5](./equations/eq5.png) 
+
+To visualize the flight trajectories predicted by this model, we integrate the spatial
+coordinates. The position of the glider on a vertical plane will be designated by
+coordinates (x, y) with respect to an inertial frame of reference, and are obtained
+from:
+
+![eq6](./equations/eq6.png)
+
+Augmenting our original two differential equations by the two equations above, we have a
+system of four first-order differential equations to solve. We will use a time-stepping
+approach, like in the previous lesson. To do so, we do need *initial values* for every
+unknown:
+
+![eq6](./equations/eq6.png)
+
+where we are now using a superscript $n$ to indicate the $n$-th value in the time
+iterations. The first differential equation, for example, gives:
+
+![eq7](./equations/eq7.png)
+
+The full system of equations discretized with Euler's method is:
+
+![eq8](./equations/eq8.png)
+
+*  Assume L/D of 5.2 (a value close to measurements in [Feng et al. 2009](./feng et al 2009-paper_airplane.pdf))
+*  For the trim velocity, v_t=5.5 m/s.
+
+**a\.** Create the ordinary differential equation function that outputs: v', theta', x', and y'
+as the function `dy=phugoid_ode(t,y)`, where `y=[v,theta,x,y]'`
+
+**b\.** Plot the position of the airplane over 20 seconds if the initial height is y(0)=2 m,
+initial position is x(0)=0 m, initial velocity is v(0)=10 m/s, and initial angle,
+theta(0)=0 rad. 
+
+  **i\.** Use an Euler approximation with timesteps of 0.1 s and 0.01 s.
+
+  **ii\.** Use ode23 
+
+  **iii\.** Include the three plots on a single graph in your README.
+
+**c\.** Group assignment for **Extra Credit**. Determine the optimal initial angle and initial speed to
+launch the paper airplane to maximize distance traveled. You can solve the path of the
+flight any way you like and you can work with up to three colleagues. 
+
+
+Submit your github repository solution to the following Google form:
+
+[https://goo.gl/forms/QpINqnnEhMceqCy02](https://goo.gl/forms/QpINqnnEhMceqCy02)
+
+Initial judging by Prof. Cooper will be based upon clarity of solution (e.g. documents, files, results in README)
+
+The top groups will be distributed to the class to be voted on
+
+**Good Luck!**
 
-### Mini-design challenge 

HW6/equations/eq1.tex

@@ -0,0 +1,4 @@
+\begin{align}
+m \frac{dv}{dt} & = - W \sin\theta - D \\
+m v \, \frac{d\theta}{dt} & = - W \cos\theta + L
+\end{align}

HW6/equations/eq2.tex

@@ -0,0 +1,4 @@
+\begin{align}
+ \frac{v'}{g} & = - \sin\theta - D/W \\
+ \frac{v}{g} \, \theta' & = - \cos\theta + L/W
+ \end{align}

HW6/equations/eq3.tex

@@ -0,0 +1 @@
+L/W=v^2/v_t^2

HW6/equations/eq4.tex

@@ -0,0 +1,4 @@
+\begin{align}
+ L &=& C_L S \times \frac{1}{2} \rho v^2 \\
+ D &=& C_D S \times \frac{1}{2} \rho v^2
+ \end{align}

HW6/equations/eq5.tex

@@ -0,0 +1,4 @@
+\begin{align}
+  v' & = - g\, \sin\theta - \frac{C_D}{C_L} \frac{g}{v_t^2} v^2 \\
+   \theta' & = - \frac{g}{v}\,\cos\theta + \frac{g}{v_t^2}\, v
+   \end{align}

HW6/equations/eq6.tex

@@ -0,0 +1,9 @@
+$$
+v(0) = v_0 \quad \text{and} \quad \theta(0) = \theta_0\\
+x(0) = x_0 \quad \text{and} \quad y(0) = y_0
+$$
+
+\begin{align}
+x'(t) & = v \cos(\theta) \\
+y'(t) & = v \sin(\theta).
+\end{align}

HW6/equations/eq7.tex

@@ -0,0 +1,2 @@
+$$\frac{v^{n+1} - v^n}{\Delta t} = - g\, \sin\theta^n - \frac{C_D}{C_L} \frac{g}{v_t^2}
+(v^n)^2$$

HW6/equations/eq8.tex

@@ -0,0 +1,8 @@
+\begin{align}
+v^{n+1} & = v^n + \Delta t \left(- g\, \sin\theta^n - \frac{C_D}{C_L} \frac{g}{v_t^2}
+(v^n)^2 \right) \\
+\theta^{n+1} & = \theta^n + \Delta t \left(- \frac{g}{v^n}\,\cos\theta^n +
+\frac{g}{v_t^2}\, v^n \right) \\
+x^{n+1} & = x^n + \Delta t \, v^n \cos\theta^n \\
+y^{n+1} & = y^n + \Delta t \, v^n \sin\theta^n.
+\end{align}