How to Use the Differential Equation Plotter -
A Brief Tutorial on Numerical Methods
What it does:
This program plots systems of up to three first order differential equations using the standard fourth order Runge-Kutta method. The independent variable is x or t and the dependent variables are x (if t is independent), y, z, u, v, and w. Initial conditions are set at t0, x0, y0, etc. The range of values plotted is tmin to tmax, xmin to xmax, ymin to ymax, etc. The step size is set at 0.001 by default but can be increased for speed or decreased for accuracy. Plotting coordinates (x vs t, y vs. x, u vs. x, etc) can also be varied.
OK fine, now what does all this mean?
How numerical methods work:
Differential euqations are solved numerically using an approximate expression for the derivative of a function to generate an iterative scheme - a method of generating new values based on previous values. The simplest method makes use of the definition of a derivative -
f'(x) = lim(f(x+h) - f(x))/h as h-->0
Taking h to be some small positive number, say 0.001, create an array of x values, x[n], starting at some value xmin and progressing up to xmax in increments of h; x = xmin, x = xmin+h, x=xmin+2h, etc. until x[final] = xmax. Now given a differential equation of the form y' = f(x, y), where f is some expression involving x and y, find an array y with values y[n] corresponding to x[n] in the array x above. Since x[n+1] = x[n] + h and y = f(x), y[n] = f(x[n]) and y[n+1] = f(x[n] + h), so if h is small enough a reasonable approximation to y' at x = x[n] is
y' = (y[n+1] - y[n])/h or
f(x,y) = (y[n+1] - y[n])/h.
Solving for y[n+1] we have:
y[n+1] = y[n] + h*f(x[n],y[n]).
Notice that given a knowledge of y[n], we can now determine y[n+1]. This is the crux of all numerical differential equation methods - given starting values for x and y (also known as initial conditions), we can now compute all subsequent values of y by literally stepping up from the initial values (hence the terms step size and initial values used in the rather cryptic "What it does" section). Obviously the closer h is to zero, the closer the approximation is to the actual derivative and the the more accurate the solution is. Unfortunately, smaller step sizes h also mean more calculations are required and the method outlined above (known as Euler's method) is notoriously inaccurate for larger step sizes, so this program uses a much more accurate algorithm known as the Runge-Kutta method. The details are too complicated to go into here but can be found in any textbook on differential equations (my personal favorite is "Elementary Differential Equations and Boundary Value Problems" by Boyce and DiPrima). Suffice it to say it is an iterative scheme in the same spirit as the Euler method outlined above. Now for a concrete example to illustrate the use of the program:
A step-by-step example - the Lorenz equations:
A study of convective motion in the atmosphere due to temperature differences led Edward Lorenz to the following third order system:
dx/dt = a(y-x),
dy/dt = rx -y -xz,
dz/dt = -bz + xy.
They may not look like much on the surface, but analysis of these equations led to an entirely new branch of mathematics - chaos theory. For the earth's atmosphere a = 10 and b = 8/3. Values for r depend on the temperature differences between different layers of the atmosphere and can assume a wide range of values; choose r = 28 for starters and plug the equations into the plotter as follows:
Enter the equations-
The Lorenz equations have t as the independent variable and x, y, and z as the dependent variables. Under "Miscellaneous Controls" choose "t, x, y, z" from "variables" by clicking on the triangle and selecting from the list that appears to the side. In the section labeled "Enter Functions:" are textboxes for x' (= dx/dt), y', z', etc. In the box labeled "x':" type in
and in the boxes labeled "y':" and "z':" type in
Determine initial values and plotting ranges-
These are usually given in the problem statement but for this problem try (x0, y0, z0) = (5, 5, 5) at t0 = 0 (enter them in the appropriate textboxes under "Initial Values"). Good plotting ranges are t: 0 to 40, x: -25 to 25, y: -25 to 25, z: -25 to 25. Enter them in the "Enter Plotting Ranges:" section. Leave the step size at 0.001 and choose "x vs. t" from "coords:". Java textboxes behave just like regular text editors so errors can be corrected using the backspace and delete keys just like in Notepad, emacs, or the DOS editor. When finished click the "Graph" button to start the plotting process (no need to hit "Enter" when typing in data; just type it in and use the mouse or "Tab" key to move from box to box). After about 5 to 10 seconds (on a 500 MHz machine) a highly chaotic graph should appear. Now click on the button on the right hand side of the "coords" box and choose "y vs. x" from the popup menu. Click "Graph" again and 40 seconds of trajectory in the x-y plane will appear. Obviously there is much more to these equations than appears on the surface! Play around with different values of r and different initial conditions and see what happens.
Higher order equations:
Equations up to third order may be handled by breaking them up into a first order system of equations. Consider the equation for the simple harmonic oscillator:
y'' + y = 0, y(0) = 3, y'(0) = 0.
Let u = y', then u' = y'' and the following system results:
y' = u, y(0) = 3
u' = -y, u(0) = 0.
Data enters in the plotter as follows:
Choose "x, y" from "variables" and click "Graph" (any unneeded fields, such as v', may be left as is). The graph should be a cosine curve of amplitude 3.
Equations as high as fourth order in x and y can be handled by using y' = u, u'=v (=y''), v'=w, and w'=f(x,y,u,v,w), where f is the expression obtained by solving for the highest order derivative and substituting u for y', v for y'', etc. By letting u=x', v=y', and w=z', systems of up to 3 second order equations may be plotted (think Langrangian dynamics!).
Things to be aware of:
The only constraint on the Runge-Kutta method (as well as most other common methods) is that you must be able to solve for the highest order derivative; i.e., put the equation in one of the forms y' = f(x, y), y'' = f(x,y,y'), or y''' = f(x,y,y', y''). Assuming this can be done there remains the problem of singular points. Consider Bessel's equation of order 0:
x2y'' + xy' + x2y = 0, y(0) = 3, y'(0) = 0.
Solving for y'' we get
y'' = -y'/x -y,
and a regular singular point can be seen at x = 0. Letting y' = u we get the system
y' = u, y(0)=3
u' = -(u/x)-y, u(0) = 0.
Assuming xmin is negative and xmax is positive, the plotter accurately depicts the solution from 0 to xmax, but in stepping back from 0 to xmin it tries to divide 0 by 0 and shows only black space. I'm currently working on an implicit method to apply in a small neighborhood about regular singular points, but time constraints may make progress slow. Please be patient and check back for updates. Any suggestions are always welcome and may be e-mailed to me at the address shown at the bottom of the calculator page.
(Note: The plotter now works around some classes of singular points, including Bessel's equation at the origin. More improvements to follow soon - I hope!)
Look for a Swing version soon.
Guy Stallings, November, 2002.