Physics 410: Solution of ODEs
Please report all errors/typos. etc to choptuik@physics.ubc.ca
Last updated November 2005
- General Reference
- Handouts:
- Part I -- Combined notes: Solution of ODEs, Orbiting
Dumbbell, Method of Lines Solution of the Wave Equation and
Shooting and The Toy Deuteron Problem (PDF).
Lecture versions
- Solution of ODEs (PDF).
- Orbiting Dumbbell (PDF).
- Method of Lines Solution of the Wave Equation (PDF).
- Shooting and the Toy Deuteron Problem (PDF).
- Part II -- Source code and related: lsoda.f through deut
(PDF)
- Source code covered in class can be found on the phys410
account on the lnx machines in sub-directories of ~phys410/ode.
- Lecture 2
- lsoda.f
ODEPACK routine for solving general systems of ODEs.
- tlsoda.f: Sample
"driver" program which uses LSODA to integrate u''(t) =
-u(t) and associated Makefile.
- chk-tlsoda.f:
Program which applies O(dt^2) finite-difference approximation of the
ODE to tlsoda results ("independent residual evaluator").
- Tlsoda: shell
script which runs tlsoda with a variety of tolerance settings,
checks one solution using "independent residual evaluation", and
demonstrates dependence of results on number of requested output times.
Output from the
script on the lnx machines.
- Plot of computed
and exact solution on x = [0 .. 10] with pure absolute error tolerance
of 1.0e-6 (gnuplot script file).
Plot of abs(error)
in computed solution using a range of tolerances (gnuplot script file).
- Sample output
illustrating the use of some utility commands: dvmesh, nf
and paste. May be useful in answering Problem 2e) of
Homework 4.
- Lectures 3 and 4
- integral: Example code illustrating use of LSODA
to compute a definite integral.
- optional-inputs.f:
Documentation of some of the optional inputs to LSODA;
information about other optional inputs can be extracted from other
parts of the comment block.
- integral.f:
Driver routine
- fcn.f: User
function which defines integrand.
- Makefile
- Sample output
on the lnx machines.
- twobody: LSODA-based code which solves
restricted 2-body gravitational problem (one body non-gravitating).
- twobody.f:
Driver routine.
- fcn.f: User
function which implements equations of motion.
- fcn.inc:
Include file which defines variables and COMMON BLOCK for communicating
additional parameters to user function.
- Makefile
- Twobody:
Shell-script which runs twobody and then produces various plots
of the output:
- Twobody 1.0 1.0d-6: (circular orbit, all LSODA
tolerances 1.0d-6)
- Test particle position (PDF).
- Total mechanical energy deviation (PDF).
- Total angular momentum deviation (PDF).
- Twobody 1.0 1.0d-10: (circular orbit, all LSODA
tolerances 1.0d-10)
- Test particle position (PDF).
- Total mechanical energy deviation (PDF).
- Total angular momentum deviation (PDF).
- Twobody 0.8 1.0d-10: (elliptical orbit, all LSODA
tolerances 1.0d-10)
- Test particle position (PDF).
- Total mechanical energy deviation (PDF).
- Total angular momentum deviation (PDF).
- Twobody-sm:
A variant of Twobody which uses supermongo
(sm)
- Twobody-sm 0.8 1.0d-10: (elliptical orbit,
all LSODA tolerances 1.0d-10)
- Test particle position (PDF).
- Total mechanical energy deviation (PDF).
- Total angular momentum deviation (PDF).
- Twobody-xvs:
A variant of Twobody which sends output directly to the xvs visualization server.
- dumb: LSODA-based code which solves orbiting dumbbell
problem
- Lecture notes (PDF)
- dumb.f: Main
program
- fcn.f: User
function which implements equations of motion.
- fcn.inc:
Include file which defines variables and COMMON BLOCK for communicating
additional parameters to user function.
- Makefile
- Dumb: Script which
runs dumb and produces plots shown below.
- Animation of typical elliptical orbit with d=0.3 (1.3MB MPEG)
- Plots of omega (angular frequency of dumbbell about its
center of mass) versus time for a circular
orbit and for an elliptical
orbit. Note the "chaotic" nature of omega in the latter case. Detail of omega for the
elliptical orbit.
- Plots illustrating energy conservation for elliptical
orbit (chaotic case). Energy
quantities computed with LSODA tolerance 1.0e-6. Top to bottom:
Translational kinetic energy, rotational kinetic energy, total energy,
gravitational potential energy. In this case, "non-conservation" of
total mechanical energy is a good sign that we need to make the error
tolerance more stringent. Energy
quantities computed with LSODA tolerance = 1.0e-12. Note that
energy conservation is improved in this case. Plot of rotational kinetic energy.
- wave: LSODA-based code which solves wave equation
using the method lines and an O(h^2) disretization of the spatial part
of the wave operator.
- Lecture notes (PDF)
- wave.f: Main
program
- fcn.f: User
function which implements equations of motion.
- fcn.inc:
Include file which defines variables and COMMON BLOCK for communicating
additional parameters to user function.
- Makefile
- Wave: Script which
runs wave and produces plot shown below. Output of script on lnx1.
- Surface plot
of wave equation solution using time symmetric initial conditions and
Dirichlet boundary conditions.
- deut: LSODA-based code which integrates
spherically-symmetric time-independent Schrodinger equation for toy
deuteron model (square well potential).
- Lecture notes (PDF)
- deut.f: Main
program
- fcn.f: User
function which implements equations of motion.
- fcn.inc:
Include file which defines variables and COMMON BLOCK for communicating
additional parameters to user function.
- Makefile
- Shoot-deut:
Script to perform bisection search ("shooting") for energy eigenvalue
using bisection search script package.
- Shoot-deut-all:
Yet another script that calls Shoot-deut to compute
wavefunctions for x0 = 2.0, 4.0, 6.0, 8.0. Note that as per the notes,
the initial bisection bracket [ELO,EHI] is to be
chosen so that u(x) is negative/positive at large x for integrations
with E=ELO/EHI respectively. It turns out that
this actually requires ELO > EHI
- Plot and detail of several trial
solutions generated via bisection search on eigenvalue, E for x0=2.
- Plot of computed
eigenfunctions for various values of parameter, x0. Note that
eigenfunctions are not normalized.
- mkplots: Script
which makes above plots using gnuplot