Physics 410: Solution of ODEs

Last updated November 13, 2000

• General Reference
  • Numerical Recipes, Chapters 16 and 17
  • ODEPACK Source Code (browse directory)
• Source code covered in class can be found on the phys410 account on sgi1 in ~phys410/ode/ex1 and ~phys410/ode/ex2
• Lectures 1 and 2 (Nov 21 and 23)
  • Handout (PS) and slides (PS).
  • lsoda.f ODEPACK routine for solving general systems of ODEs.
  • tlsoda.f Sample driver routine which uses LSODA to integrate u''(x) = -u(x). Makefile and sample usage on SGIs. Plot of computed and exact solution on x = [0 .. 15] with all error tolerances set to 1.0d-6. Plot of error in computed solution. Note that this example also requires the source file dvfrom.f.
• Lectures 2, 3 and 4 (Nov 23 and 28)
  • Handout (PS) and slides (PS).
  • integral: Example code illustrating use of LSODA to compute a definite integral.
    • integral.f: Driver routine
    • fcn.f: User function which defines integrand.
    • Makefile
    • Sample output on SGIs.
  • 2body: LSODA-based code which solves restricted 2-body gravitational problem (one body non-gravitating).
    • 2body.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
    • Doit: Shell-script "front-end" to 2body which illustrates many shell programming techniques. Sample output:
      • Doit 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).
        • Output from script with tracing enabled (-x flag).
      • Doit 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).
      • Doit 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).
    • Doit-s: The same script, but with the majority of the comments removed.
    • Doit-sm: A variant of Doit which uses supermongo (sm) to do the plotting. Sample output
      • Doit-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).
    • Doit-jser: A variant of Doit which sends output directly to the scivis/jser visualization server using the jv1 client program.
• Lectures 4 and 5 (Nov 28)
  • Handout (Part 1 PS and Part 2 PS) and slides (Part 1 PS and Part 2 PS)
  • dumb: LSODA-based code which solves orbiting dumbbell problem
    • dumb.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
    • 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.