Physics 410 - Computational Physics - UBC - Fall 2002

Physics 410: Solution of ODEs

Please report all errors/typos. etc to [email protected]

Last updated November 2002

General Reference
- Numerical Recipes, Chapters 16 and 17
- ODEPACK Source Code (browse directory)
Handouts:
- Part I -- lsoda.f through deut (PS)
- Part II -- Combined notes: Orbiting Dumbbell, Method of Lines Solution of the Wave Equation and Shooting and The Toy Deuteron Problem (PS)
Source code covered in class can be found on the phys410 account on the lnx machines in sub-directories of ~phys410/ode.
Lecture 3 (Nov 14)
- 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.
Lectures 3, 4 and 5 (Nov 14, 19 and 21)
- 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 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
  - 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.
  - 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.
  - 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).
  - 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.
  - 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
Supplementary Material for Homework 5 (source code available via ~phys410/hw5)
- gpwave.f: This program illustrates the output of time series of one-dimensional profiles (i.e. f(x,t)) for subsequent plotting as a surface plot using gnuplot. Sample usage on lnx1; gnuplot script file file, and surface plot (PS). Use help splot in gnuplot for more information on making surface plots.