Physics 410: Solution of ODEs

Last updated November 2004

• General Reference
  • Numerical Recipes, Chapters 16 and 17
  • ODEPACK Source Code (browse directory)
• 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 (PS).
    Lecture versions
    1. Solution of ODEs (PS).
    2. Orbiting Dumbbell (PS).
    3. Method of Lines Solution of the Wave Equation (PS).
    4. Shooting and the Toy Deuteron Problem (PS).
  • Part II -- Source code and related: lsoda.f through deut (PS)
• Source code covered in class can be found on the phys410 account on the lnx machines in sub-directories of ~phys410/ode.
• Lecture 2 (Nov 23)
  • 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 2, 3 and 4 (Nov 23, 25 and 30)
  • 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 (PS).
        • Total mechanical energy deviation (PS).
        • Total angular momentum deviation (PS).
      • Twobody 1.0 1.0d-10: (circular orbit, all LSODA tolerances 1.0d-10)
        • Test particle position (PS).
        • Total mechanical energy deviation (PS).
        • Total angular momentum deviation (PS).
      • Twobody 0.8 1.0d-10: (elliptical orbit, all LSODA tolerances 1.0d-10)
        • Test particle position (PS).
        • Total mechanical energy deviation (PS).
        • Total angular momentum deviation (PS).
    • 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 (PS).
        • Total mechanical energy deviation (PS).
        • Total angular momentum deviation (PS).
    • 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 (PS)
    • 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 (PS)
    • 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 (PS)
    • 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 [E_LO,E_HI] is to be chosen so that u(x) is negative/positive at large x for integrations with E=E_LO/E_HI respectively. It turns out that this actually requires E_LO > E_HI
    • 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