Many physical problems are modeled with differential equations, and today we'll consider a simple one-dimensional problem that is the starting point for understanding the higher dimensional cases that arise in practice. 1D problems are often used to study higher dimensional problems, because they are simpler and less expensive to solve numerically.
A variety of methods are used to solve differential equations, and here we'll look at an iterative method used to solve the boundary value problem given by the following ordinary differential equation:
We will solve this equation on the unit interval [0,1] subject to boundary conditions
We will use the method of finite differences to solve the problem; we discretize the solution u by considering the values of u only at discrete points spaced regularly along the unit interval.
We divide the unit interval into N-1 subintervals using N points positioned uniformly along the line at positions x=i×h, i=0, 1, 2, ..., N-1. Let u[i] represent the value of the approximate solution at point x = i×h, u(i×h), where h = 1/(N-1) and is called the mesh spacing
i=0 i=1 i=2 i=N-1
|-----|-----|------| ..... |-----|-----|
^ ^
| |
x=0 x=h x=2h x=1-h x=1
We use a Taylor series to approximate the second derivative u''(x) on the left hand side of equation (ODE) as
u''(i*h) ~ ( u[i-1] - 2u[i] + u[i+1] )/h2 (discrete 2nd derivative)
Setting this equation equal to the right hand side of equation (ODE), and rearranging terms, we come up with the following iterative algorithm, known as Gauss-Seidel method.
repeat {
for i = 1 to N-2
u[i]= (u[i-1] + u[i+1] - 8*h*h ) / 2
}
until ( Error(u,uExact) < ε);
where uExact[ ] contains the exact solution.
The method updates each value of u[i] as a function of its nearest neighbors:
u[i-1] u[i] u[i+1] |-----------|-----------| x = (i-1)*h i*h (i+1)*h
It can be shown that the Gauss-Seidel method will converge an approximate solution to u. The smaller we make h, the more accurate our answer will be, but at a greater cost, since N increases as h-1. If you are interested in the details, consult Numerical Recipes in C, by Press et al., Cambridge University Press, which is also available on line (see Chapter 19, Relaxation Methods for Boundary Value Problems, pp 863-4), or a good book on numerical analysis.
The values of u[i] when i is 0 or N-1 are boundary data, and are given as input; these are defined to be 0 in our case. The evolving solution is represented on the interior N-2-element sub-array, u[1] through u[N-2]. Beginning with an initial guess we update the solution according to the right hand side of the equation (ODE) until we are satisfied with the answer. We are finished computing when the error drops below the user-specified value ε. It can be shown that the solution will converge in O(N2) iterations. The more accurate the answer required, the greater the cost. The cost of computing the solution increases with N and with decreasing tolerance.
In general we do not know the exact solution so we must estimate the error to determine when our computed solution has reached (or converged to) the desired degree of accuracy (Our simple linear ODE has a known exact solution: u(x) = 4x(x-1).). The residual is a useful way to estimate the error. To compute the residual, we plug in our computed solution into the original equation, and note the difference between the resultant left and right hand sides, we do this for all x, where D2(xi) is the discrete approximation to the second derivative of the computed solution u at x=xi, and is given the formula displayed above:
r(xi) = D2(xi) - 8
We treat the difference as a vector, one for each value of x. The residual is the length of this vector normalized by the number of elements in the solution: √∑r[i]2 / n
One last detail remains. What values should we choose for the initial values of the solution array u[ ]? We usually come up with some reasonable guess to the solution; the better the guess the more quickly the answer will converge, and fewer iterations will be required. Coming up with a good guess depends on the problem. In our case we will choose the guess to be a solution of all ones, except for the boundary points which will be set to zero as described above.
The trouble with the above numerical algorithm is that it doesn't parallelize. Each value of u[i] depends on all values of u[ ] to the left and to the right:
for i in 1 : N-2
u[i]= (u[i-1] + u[i+1] - 8*h*h) / 2
end for
One way to eliminate the problem is to rename the left hand side of the update formula. We maintain both and old and new meshes, storing updates into the new mesh. This results in an algorithm known as Jacobi's method.
It can be shown the the Jacobi method will eventually reach
(converge to) the exact solution:
repeat s
for i in 1:N-2
ui(s) = (ui-1(s-1) + ui+1(s-1) )/ 4
end for
until |u(s) - u(s-1)| is small enough
where ui(s-1) is the value of the mesh at (i) at iteration s-1. Each value in the solution is the average of the 2 nearest neighbors in the "old" mesh. To save space, we swap old and new meshes after each iteration.
Another approach is to divide the unknowns into two groups---red and black (sometimes called odd and even) then we can decouple the points into disjoint sets, and update them in any desired order.
@ -> boundary @----o----o----o----o----o----o----o----o----o----o----@
The algorithm is as follows:
repeat
// Update red points
for i in 1:N-2, i is even
u[i] = (u[i-1] + u[i+1] - 8*h*h ) / 2
endfor
// Update black points
for i in 1:N-2, i is odd
u[i] = (u[i-1] + u[i+1] - 8*h*h ) / 2
endfor
until error is small enough.
With Red/Black ordering the solution may converge more quickly than with the usual ordering. However, the computation within the inner loop is a bit more expensive since it must access (read) each point in the solution twice per sweep (Why is this so?). However, the benefit of this algorithm is that it parallelizes. This technique valuable in removing data dependencies in a variety of algorithms.
Although Red/black ordering offers improved performance over point Jacobi, it doesn't improve on the slow asymptotic running time: O(N2) iterations for an N=n×n mesh!
The problem is that the lower frequency components of the error converge very slowly and so we need to let these errors diminish over a large number of iterations. However, the high frequency components converge very quickly, often in just a few iterations. Herein lies a strategy for reducing the running time to a constant number of iterations, and algorithm that uses this strategy known as Multigrid.