Index of /www/kankel/ph567/examples/Octave/ODE

	Name	Last modified	Size	Description
	Parent Directory		-
	233px-Messier_53_HST.jpg	2022-03-28 11:25	30K
	DoublePendulum/	2022-04-04 21:32	-
	TwoStep.png	2022-03-23 13:57	117K
	fft_greens_function/	2011-03-30 12:59	-
	hierarchy.fig	2020-02-14 13:18	4.1K
	hierarchy.fig.bak	2020-02-14 13:18	4.0K
	hierarchy.pdf	2014-02-24 09:36	20K
	lorenz/	2011-02-17 21:30	-
	plots/	2009-03-09 16:26	-
	programs/	2018-03-10 02:46	-
	shuffling/	2020-02-14 13:18	-

Ordinary Differential Equations

Introduction

Outline:

1. Fixed Step Size
1. Euler's method
2. Runge-Kutta methods
2. Quadrature
3. Step doubling
4. Variable step size
5. Advanced solvers (LSODE)
6. Sensitivity to initial conditions; Chaos
7. Solution of stiff ODEs (and numerical differentiation)
8. Shooting \& root finding
9. Green's function methods in the Fourier domain

Along the way, as an example, I will be developing my own toolkit for solving ODEs. This is not intended to replace the high quality ODE solvers available in Octave or Matlab. Instead, my purpose is to demonstrate principles of modular code design and good documentation. I will also demonstrate the use of recursion.

We will begin with techniques using a fixed step size. Before looking at the examples in this section, read the lecture notes.

Question: In physics, an ansatz (German 'answer') is an inspired guess that is verified afterward. What motivation can you find for the ansatz proposed on page 4?

Fixed Step Size

Euler's Method

Source code:

• euler.m
• ode_fixed.m
• harmonic.m
• euler_test1.m

Plots:

• Harmonic oscillator y(t)
• Phase space trajectory

RK2

Source code:

• rk2.m
• rk2_test1.m
• goofy.m
• rk2_test2.m

Plots:

• Harmonic oscillator y(t)
• Phase space trajectory
• Nonlinear oscillator

Quadrature

Step Doubling

First, visit the lecture notes.

Source code:

• twostep.m
• Apply to SHO: twostep_test1.m
• Error estimation: twostep_test2.m

Plots:

• Harmonic oscillator y(t)
• Phase space trajectory
• Error estimation

Variable Step Size

Having demonstrated that the step double works, we now build a fully functional adjustable spanner adjustable step size solver.

Source code:

• The engine that makes it work: adaptive.m
• The front end: adaptive_stepper.m
• Apply to SHO: adaptive_test1.m

Plots:

• SHO phase space
• SHO convergence

Hi. My name is Isaac, and I have a problem.

The N-body problem is notoriously difficult. Let's try it! But first, read the lecture notes.

Source code:

• adaptive_stepper.m
• System of ODEs: Nbody.m
• NbodyInvariants.m
• Nbody_test1.m
• Sensitivity to initial conditions: Nbody_test2.m
• Comparison to LSODE: Nbody_test3.m

Plots:

• Test 1
• Test 2: Sensitivity to initial conditions
• Test 3: LSODE Comparison

My solver works surprisingly well, but it's not a professional grade tool. It is inefficient for many reasons (can you name a few?). LSODE is more efficient and more robust, but even LSODE quickly reaches its limits. For this problem, the best tools are symplectic integrators. Ask Neil!

It might be wortwhile now to review the hierarchy of modules developed thus far.

Chaos

Chaos is mainly about sensitivity to initial conditions. That is, very slightly different initial conditions can give rise to qualitatively different solutions. See the lecture notes.

Source code:

• ODEs: lorenz.m
• Demonstration: sensitivity.m

Plots:

• Lorenz attractor plot
• Lorenz attractor movie
• Sensitivity to initial conditions

Shuffling cards: an analogy to chaos

Conceptually, we can view a chaotic system as 'shuffling' in phase space. Here we will explore this analogy with a simple iterative map. The shuffling idea can help us to place a lower bound on the phase space dimensionality required to accommodate chaotic behavior.

• Conceptual sketch
• Source code: iterative_map.m
• Bifurcation diagram

Quantum Chaos

What happens when you add a quantum-scale perturbation to a classical, Hamiltonian, macroscopic chaotic system?

While we are on the topic quantum chaos, we may optionally investigate dynamical billiards. The quantum version of this system can be solved using the PWDM technique. Unfortunately, this has nothing to do with ODEs.

Stiff Equations

If the equations of motion turn out to be stiff, your ODE solver might be DOA! To solve this problem, you'll need a time machine.

• Solver: backward_euler.m
• Improved front end: adaptive_stepper2.m
• EOM: stiffeq.m
• Jacobian: stiffeq_jacobian.m
• Demonstration: stiff_test1.m
• Results!

Here again is my hierarchy of modules. You now know when & why to use the stiff solvers.