where is the independent variable, are the dependent variable, and underlines denote vectors.

This looks like a first-order ODE, but it can represent a higher order ODE with a simple "chain" of variables. E.g., can be implemented by defining , and

We can also eliminate explicit dependence by making it one of the elements of .

Note: I'm going to drop the underlines now, but and are still vectors.

- We want to find approximate numeric solutions of some accuracy without already knowing the analytic solution.
- For initial value problems, we do this by integrating the ODE in discrete steps in . The step size is conventionally called .
- The easiest starting point is the Taylor series.

Euler 1st order:

`***` If we require that only the function be specified
in analytic form, not its derivatives, then has to be estimated from
previous or future values of in some clever way.

where is the Euler's method approximation for .

This is equivalent to

where

General principles:

- Don't embed an algorithm everywhere you need it. That forces you to rewrite it repeatedly.
- Write once, reuse often!
- This is a little more work, but usually pays for itself quickly.
- The little bit more work: you have to design an API for your algorithm.

The API can be very simple: just one function to advance one step.

The ODE step function takes three arguments:

- The starting value of . (A
`vector<double>`.) - The function for . (A pointer to a function.)
- The size of the step .

There are two options for how to return the updated value:

- Actually return the new value.
- Update the existing "in place" (pass by reference).

The function takes a little more thought.

The straight-forward way: pass-by-value, return a value:

vector <double> f( vector<double> y );

Pass-by-reference (in and out):

void f(const vector<double> &y, vector<double> &f_out);

The second way is up to a factor of 2 faster because it requires less copying of data between memory locations.

typedef void f_vect_vect_t(const vector<double> &y, vector<double> &f_out); void rk4step(vector<double> &y, f_vect_vect_t * fptr, double h);

Good "namespace" citizenship:

- Choose names unlikely to be used elsewhere. (E.g., "f" would be a bad name for a globally-defined function or type.)
- C++
`namespaces`were designed to help with this.

- Write test cases for your functions.
- You need a test independent from your main application!
- Don't be fooled into thinking you can always recognize a wrong answer when you see it. If you already knew the answer to your main problem, you wouldn't be writing this program in the first place.

We know the solutions to and to . Let's test the Euler ODE solver with that case.

Code: see course web page.

Position vs time:

It's spiralling outward. Energy is not being conserved!

Momentarily using the underline notation again, the vector
is a vector in *phase space*.
It is supposed to follow a closed path, but the finite steps in the
1-st order Euler ODE solver don't allow that.

The harmonic oscillator and orbit problems, and many others in science in engineering, are examples of solutions of Hamiltonian equations:

The solutions should conserve . The Euler and Runge-Kutta schemes don't.

The simplest (and very common) cases arise from , so that and .

A 1st-order symplectic version of Euler's equation:

A beautiful 2nd-order symplectic method:

Note these symplectic methods need to distinguish between and .

- Implement the 2-nd order symplectic integrator and the 4-th order Runge-Kutta integrator.
- Test them with oscillator and orbit equations.
- You can reuse as much of the code on the course web page as you like.