<- previous index next ->
Quick and dirty numerical solution
Using the notation y for y(x), y' for d y(x)/dx, y'' for d^2 y(x)/dx^2 etc.
To solve a(x) y''' + b(x) y'' + c(x) y' + d(x) y + e(x) = 0
for y at one or more values of x,
with initial condition constants c1, c2, c3 at x0:
y(x0) = c1 y'(x0) = c2 y''(x0) = c3
Create a function that computes y'''
f(ypp, yp, y, x) = -(b(x)*ypp + c(x)*yp + d(x)*y + e(x))/a(x)
a(x), b(x), c(x), d(x) and e(x) may be any functions of x and
that includes constants and zero. a(x) must not equal zero in
the range of x0 .. desired largest x.
Initialize x = x0
y = c1
yp = c2
ypp = c3
loop: use some method with the function, f, to compute the next ypp
yppp=f(ypp, yp, y, x)
ypp = ypp + h*yppp first order or use some higher order
yp = yp + h*ypp or some higher order method
y = y + h*yp or some higher order method
x = x +h increment x and loop
Quit when you have solved for y at the desired x
You may save previous values of x and y, interpolate at desired
value of x to get y.
You may run again with h half as big, keep halving h until you get
approximately the same result.
Quick and dirty finished! This actually works sometimes.
Better numerical solutions
To see some very accurate solutions to very simple ODE's
Solve y' = y y(0)=1 which is just y(x)= exp(x)
ode_exp.c
ode_exp_c.out note very small h is worse
To see some very accurate solutions to simple second order ODE's
Solve y'' = -y y(0)=1 y'(0)=0 which is just y(x)= cos(x)
ode_cos.c
ode_cos_c.out many cycles
An example of an unstable ODE is Hopf:
hopf_ode.c
A three body problem, using simple integration, shows the
sling-shot effect of gravity when bodies get close.
This is numerically very unstable.
body3.c needs OpenGL to compile and link.
Hopefully, it can be demonstrated running in class.
body3.java plain Java, different version.
Hopefully, it can be demonstrated running in class.
For more typical ordinary differential equations there are some
classic methods. A common higher order method is Runge-Kutta fourth
order.
Given y' = f(x,y) and the initial values of x and y
L: k1 = h * f(x,y)
k2 = h * f(x+h/2.0,y+k1/2.0)
k3 = h * f(x+h/2.0,y+k2/2.0)
k4 = h * f(x+h,y+k3)
y = y + (k1+2.0*k2+2.0*k3+k4)/6.0
x = x + h
quit when x>your_biggest_x
loop to L:
When the solution y(x) is a fourth order polynomial, or lower
order polynomial, the solution will be computed with no truncation
error, yet may have some roundoff error. Very large step sizes, h,
may be used for this very special case.
In general, the solution will not be accurately approximated by
a low order polynomial. Thus, even the Runge-Kutta method may
require a very small step size in order to compute an accurate
solution. Because very small step sizes result in long computations
and may have large accumulations of roundoff error, variable
step size methods are often used.
Typically a maximum step size and a minimum step size are set.
Starting with some intermediate step size, h, the computation
proceeds as shown below.
Given y' = f(x,y) and the initial values of x and y
L: y1 = y + dy1 using some method with step size h compute dy1
the value of y1 is at x = x + h
y2a = y + dy2a using same method with step size h/2 compute dy2a
the value of y2a is at x = x + h/2
y2 = y2a + dy2 using same method, from y2a at x = x + h/2
with step size h/2 compute dy2
the value of y2 is at x = x + h
abs(y1-y2)>tolerance h = h/2, loop to L:
abs(y1-y2) < tolerance/4 h = 2 * h
y = y2
x = x + h
loop to L: until x > final x
There are many variations on the variable step size method.
The above is one simple version, partially demonstrated by:
FitzHugh-Nagumo equations system of two ODE's
v' = v - v^3/3.0 + r initial v = -1.0
r' = -(v - 0.2 - 0.2*r) initial r = 1.0
initial h = 0.001
t in 0 .. 25
fitz.c just decreasing step size
fitz.out about equally spaced output
fitz.m easier to understand, MatLab
Of course, we can let MatLab solve the same ODE system
fitz2.m using MatLab ode45
The ode45 in MatLab uses a fourth order Runge-Kutta and then
a fifth order Runge-Kutta and compares the result. A neat formulation
is used to minimize the number of required function evaluations.
The fifth order method reuses some of the fourth order evaluations.
Written in "C" from p355 of textbook
fitz2.c similar to fitz.c
fitz2.out close to same results
Another coding of fourth order Runge-Kutta known as Gill's method
is the (subroutine, method, function) Runge given in a few languages.
The user provided code comes first, then Runge.
This code solves a system of first order differential equations and
thus can some higher order differential equations also, as shown above.
The system of equations is N equations in N dependent variables, Y,
with one independent variable X.
Y_i = Y'_i(X, Y_1, Y_2, ... , Y_N) for i=1,N
The user provides the code for the functions Y'_i and places
the results in the F array. X goes from initial value, set by user,
to XLIM final value set by user. The user initializes the Y_i.
sim_difeq.f90 Fortran version
sim_difeq_f90.out
sim_difeq.java Java version
sim_difeq_java.out close to same results
sim_difeq.c C version
sim_difeq_c.out close to same results
sim_difeq.m MatLab version (main)
ode_test1.m funct computes derivatives
ode_test4.m funct computes derivatives
<- previous index next ->