Day 28, Mat360: Numerical Analysis

Today:

• Announcements:
  • You have a homework due today, over section 4.4.

    This is your last "homework" to turn in; any other assignments will be part of your final takehome (thinking of 6.2, 6.3). And these are given on your assignment page.

  • Exam 2 corrections returned

  • Still haven't inspected your sigmacs. Sorry about that. (By the way, did you remember to attach your blank check, and write your mother's maiden name? :)

  • Next time we've got a review for the final; then Wednesday is project day.

  • Reminder that Danny is speaking Friday, 3:00pm, in MP461.

• Let's recall what we did last time:
  • Learned the general method for solving differential equations....:)
  • We had a quick overview of ODEs, and the strategies we take to approximate their solutions:
    • If you're lucky, you can use the general method.
    • Replace the offending ODE by a simpler one.
    • Use numerical methods on the actual ODE.

  • Then we discussed the Taylor methods for approximating solutions, starting with
    • Tangent lines (Euler, or T1)
    • Tangent quadratics (T2)
    • Tangent cubics (T3)
    The higher order methods (T2, T3) involve computing derivatives, which we are loathe to do.

• Today: we end on Runge-Kutta methods
  • Before we get to those however, let's talk about writing second order ODEs as first order systems (p. 212). Let's look at some examples, revisiting the Taylor examples from last time:
    • Mathematica code for the example of page 209 -- and more! (pdf).

  • In section 6.3, our author takes a clever approach to create these "foundational" rk-methods, using what we've learned in chapter 4 on integration to create some methods which (approach) the classic numerical scheme, Runge-Kutta-4.

    The idea is to approach solving the ODE by using integration, which he observes tends to be a very stable operation.

    An ODE of the form $\dot{y}=f(t,y)$ with $y(t_0)=y_0$ has an exact solution that can be written in terms of an integral: $$\begin{array}{}& {\int }_{t_0}^{\tilde{t}}\dot{y}dt={\int }_{t_0}^{\tilde{t}}f(t,y)dt⟶y(\tilde{t})-y(t_0)={\int }_{t_0}^{\tilde{t}}f(t,y)dt⟶y(\tilde{t})=y(t_0)+{\int }_{t_0}^{\tilde{t}}f(t,y)dt.\end{array}$$

    Trouble is, we've got no idea what that $y$ is in the integrand; in fact, that's what we want!:)

  • We turn this into a method once we figure out how to take a small step, $$y(t_1)=y(t_0)+{\int }_{t_0}^{t_1}f(t,y)dt$$ and then do it again and again: $$y(t_{i+1})=y(t_i)+{\int }_{t_i}^{t_{i+1}}f(t,y)dt.$$

  • The next step is to integrate using our various schemes:
    • "euler" -- note its stencil (p. 215): what integration method is this?
    • trapezoidal
    • midpoint (our author doesn't include this one in the discussion, but leaves it for the exercises).
    • simpson
    • Two additional stencil-based schemes.
    • One interesting question is whether it would be easy to adapt our atraptive method to this process, and so create an adaptive ODE solver.

  • The "Euler" method derived from the stencil in the usual way (undetermined coefficients -- let's review that on p. 215) leads to $$y_{i+1}=y_i+hf(t_i,y_i)$$

    Easy peasy!

  • Beyond the "Euler" method, however, these methods lead to an implicit equation for $y(t_{i+1})$, e.g. $$y_{i+1}=y_i+\frac{h}{2}[f(t_i,y_i)+f(t_{i+1},y_{i+1})].$$ for the trapezoidal method, so we need to use some sort of trick...

    $$y_{i+1}=y_i+\frac{h}{2}[f(t_i,y_i)+f(t_{i+1},y_i+h\cdot f(t_i,y_i))].$$

  • And for the "Simpson" method we really struggle to figure out how to approximate those intervening values of $y$: $$y_{i+1}=y_i+\frac{h}{6}[f(t_i,y_i)+4f(t_{i+1/2},y_{i+1/2})+f(t_{i+1},y_{i+1})].$$

    $$\begin{array}{rcl}{k}_1& =& f(t_i,y_i)\\ k_2& =& f(t_i+\frac{h}{2},y_i+\frac{h}{2}{k}_1)\\ k_3& =& f(t_{i+1},y_i+h{k}_2)\\ y_{i+1}& =& y_i+\frac{h}{6}[k_1+4k_2+k_3].\end{array}$$

    The choice above is suboptimal as the author emphasizes in this section -- there's a failure of the method to "live up to its hype".

  • RK4 is similar, and derived in the next section on page 231 in a Crumpet: it splits the middle term, and approximates each in a different way! $$y_{i+1}=y_i+\frac{h}{6}[f(t_i,y_i)+2f(t_{i+1/2},y_{i+1/2})+2f(t_{i+1/2},y_{i+1/2})+f(t_{i+1},y_{i+1})].$$

    $$\begin{array}{rcl}{k}_1& =& f(t_i,y_i)\\ k_2& =& f(t_i+\frac{h}{2},y_i+\frac{h}{2}{k}_1)\\ k_3& =& f(t_i+\frac{h}{2},y_i+\frac{h}{2}{k}_2)\\ k_4& =& f(t_{i+1},y_i+h{k}_3)\\ y_{i+1}& =& y_i+\frac{h}{6}[k_1+2k_2+2k_3+k_4].\end{array}$$

    These choices are motivated by matching terms in Taylor series expansions, as one can see in the Crumpet.

  • Let's look at some examples of the use of these methods (including RK-4):
    • Mathematica code for the example tables (pdf).
    • Exercises #1-3ac

• Links:
  • Our Textbook: Tea Time Numerical Analysis (Leon Q. Brin)

  • Some old lisp routines illustrating the connection between adaptive quadrature and adaptive plotting:
    • Here's my lisp version of the plotting routine, then turned into an adaptive quadrature scheme:
    • Here's a similar lisp version of a local adaptive quadrature routine, which is then turned into an adaptive plotting scheme:

      It's not quite what you're asked to do in this algorithm, however, since you're to keep track of the errors, and keep poking at the largest:

    • Here's the global adaptive quadrature routine, which is then turned into an adaptive plotting scheme:

  • A paper my dad and I wrote about integration techniques: including lots of details!:)
  • Some notes
  • An interesting alternative derivation technique. Here, for example, is a four-point forward difference formula for the derivative, obtained using the Newton polynomial interpolator (notice, however, that the $x$ values in the numerator should actually be $y$ values -- function values):
    

    But the idea is really swell, and I am delighted by this method!

  • An example of cubic splines: how do we intepolate two points and their slopes using a cubic? We can derive the so-called "Hermite spline" as a variant of Lagrange interpolation, featuring "get out of the way" functions!

    Here's another picture of how we might build one of these.

    More on my implementation of this in Mathematica.

  • Homework solutions:
    • 1.1 (Mathematica)
    • 1.2 (Mathematica)
    • 1.3 (Mathematica)
    • 2.1 (Mathematica)
    • 2.2 (Mathematica)
    • 2.3 (Mathematica)
    • 2.4 (Mathematica)
    • 3.2 (Mathematica)
    • 3.3 (Mathematica)
    • 4.1 (Mathematica)

  • My Mathematica code for Section 2.6. Muller code, too.
  • My Octave Project for 2.6 homework.
  • Octave project for 2.5
  • Divided Differences. I use these in my Muller code, and they will prove useful when it comes to interpolating polynomials.
  • Author's octave code for these complex fractal convergence diagrams
  • Newton's method in Mathematica
  • fixed point with Aitken's acceleration in Mathematica
  • Some nice notes (and code) snippets on fixed point iteration, including these acceleration techniques.
  • fixed point in Mathematica
  • Check out how my bisection works, and how it is applied to arcsine....
  • Various sins....
    • Sins in Mathematica
    • Sins in Octave
    • ArcSins in Mathematica