<- previous index next ->
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
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
the first order method is yppp=f(ypp, yp, y, x)
ypp = ypp + h*yppp
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 y at the desired x
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.
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
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 minimum step size and a maximum 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
<- previous index next ->