<- previous index next ->
We now extend the analysis and solution to three dimensions and
unequal step sizes.
We are told there is an unknown function u(x,y,z) that satisfies
the partial differential equation:
(read 'd' as the partial derivative symbol, '^' as exponent)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/dz^2 = c(x,y,z)
and we are given the function c(x,y,z)
and we are given the values of u(x,y,z) on the boundary
of the domain of x,y,z in xmin,ymin,zmin to xmax,ymax,zmax
To compute the numerical solution of u(x,y,z) at a set of points
we choose the number of points in the x,y and z directions, nx,ny,nz.
Thus we have a grid of cells typically programmed as a
three dimensional array.
From the above we get a step size hx = (xmax-xmin)/(nx-1) ,
hy = (ymax-ymin)/(ny-1) and hz = (zmax-zmin)/(nz-1).
If we need u(x,y,z) at some point not on our solution grid,
we can use three dimensional interpolation to find the desired value.
Now, write the differential equation as a second order
numerical difference equation.
(Numerical derivatives from nderiv.out order 2, points 3, term 1)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/d^2z =
1/h^2(u(x-h,y,z)-2u(x,y,z)+u(x+h,y,z)) +
1/h^2(u(x,y-h,z)-2u(x,y,z)+u(x,y+h,z)) +
1/h^2(u(x,y,z-h)-2u(x,y,z)+u(x,y,z+h)) + O(h^2)
Setting the above equal to c(x,y,z)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/d^2z = c(x,y,z) becomes
1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +
1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +
1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) = c(x,y,z)
now solve for new u(x,y,z), based on neighboring grid points
u(x,y,z) = ((u(x-hx,y,z)+u(x+hx,y,z))/hx^2 +
(u(x,y-hy,z)+u(x,y+hy,z))/hy^2 +
(u(x,y,z-hz)+u(x,y,z+hz))/hz^2 - c(x,y,z)) /
(2/hx^2 + 2/hy^2 + 2/hz^2)
The formula for u(x,y,z) is programmed based on indexes i,j,k
such that x=xmin+i*hx, y=ymin+j*hy, z=zmin+k*hz
for i=1,nx-2 j=1,ny-2 k=1,nz-2
Note that the boundary values are set and unchanged for the six
sides i=0 for j=0,ny-1 k=0,nz-1 i=nx-1 for j=0,ny-1 k=0,nz-1
j=0 for i=0,nx-1 k=0,nz-1 j=ny-1 for i=0,nx-1 k=0,nz-1
k=0 for i=0,nx-1 j=0,nx-1 k=nz-1 for i=0,nx-1 j=0,ny-1
and the interior cells are set to some value, typically zero.
Then iterate, computing new interior cells based on the previous
set of cells including the boundary.
Compute the absolute value of the difference between each new
interior cell and the previous interior cell. A stable problem
will result in this value getting monotonically smaller.
The stopping condition may be based on the above difference
value or on the number of iterations.
To demonstrate the above development, a known u(x,y,z) is used.
The boundary values are set and the iterations computed.
Since the solution is known, the error can be computed at
each cell and reported for each iteration.
At several iterations, the computed solution and the errors
from the known exact solution are printed.
The same code has been programmed in "C", Java, Fortran 95 and Ada 95
as shown below with file types .c, .java, .f90 and .adb:
"C" source code pde3b.c
output of pde3b.c
Fortran 95 source code pde3b.f90
output of pde3b.f90
Java source code pde3b.java
output of pde3b.java
Ada 95 source code pde3b.adb
output of pde3b.adb
Matlab 7 source code pde3b.m
output of pde3b.m
It should be noted that a direct solution is possible.
The iterative solution comes from a large set of simultaneous
equations that may be able to be solved accurately. For the
same three dimensional problem above, observe the development
of the set of simultaneous equations:
We are told there is an unknown function u(x,y,z) that satisfies
the partial differential equation:
(read 'd' as the partial derivative symbol, '^' as exponent)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/dz^2 = c(x,y,z)
and we are given the function c(x,y,z)
and we are given the values of u(x,y,z) on the boundary ub(x,y,z)
of the domain of x,y,z in xmin,ymin,zmin to xmax,ymax,zmax
To compute the numerical solution of u(x,y,z) at a set of points
we choose the number of points in the x,y and z directions, nx,ny,nz.
Thus we have a grid of cells typically programmed as a
three dimensional array.
From the above we get a step size hx = (xmax-ymin)/(nx-1) ,
hy = (ymax-ymin)/(ny-1) and hz = (xmax-zmin)/(nz-1).
If we need u(x,y,z) at some point not on our solution grid,
we can use three dimensional interpolation to find the desired value.
Now, write the differential equation as a second order
numerical difference equation.
(Numerical derivatives from nderiv.out order 2, points 3, term 1)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/d^2z =
1/h^2(u(x-h,y,z)-2u(x,y,z)+u(x+h,y,z)) +
1/h^2(u(x,y-h,z)-2u(x,y,z)+u(x,y+h,z)) +
1/h^2(u(x,y,z-h)-2u(x,y,z)+u(x,y,z+h)) + O(h^2)
Setting the above equal to c(x,y,z)
d^2u(x,y,z)/dx^2 + d^2u(x,y,z)/dy^2 + d^2u(x,y,z)/d^2z = c(x,y,z) becomes
1/hx^2(u(x-hx,y,z)+u(x+hx,y,z)-2u(x,y,z)) +
1/hy^2(u(x,y-hy,z)+u(x,y+hy,z)-2u(x,y,z)) +
1/hz^2(u(x,y,z-hz)+u(x,y,z+hz)-2u(x,y,z)) = c(x,y,z)
Now, rearrange terms to collect the u(x,y,z) coefficients:
u(x-hx,y,z)/hx^2 + u(x+hx,y,z)/hx^2 +
u(x,y-hy,z)/hy^2 + u(x,y+hy,z)/hy^2 +
u(x,y,z-hz)/hz^2 + u(x,y,z+hz)/hz^2 -
u(x,y,z)*(2/hx^2 + 2/hy^2 + 2/hz^2) = c(x,y,z)
We build a linear system of simultaneous equations M U = V
by computing the values of the matrix M and right hand side vector V.
(The system of equations is linear if the differential equation is linear,
otherwise a system of nonlinear equations must be solved.)
The linear subscripts of U corresponding to u(i,j,k) is
m = (i-1)*(ny-2)*(nz-2) + (j-1)*(nz-2) + k
and the rows and columns and V use the same linear subscripts.
Now, using indexes rather than x, x+hx, x-hx,
solve system of linear equations for u(i,j,k).
u(i,j,k)=u(xmin+i*hx, ymin+j*hz, zmin+k*hz)
u(i-1,j,k)/hx^2 + u(i+1,j,k)/hx^2 +
u(i,j-1,k)/hy^2 + u(i,j+1,k)/hy^2 +
u(i,j,k-1)/hz^2 + u(i,j,k+1)/hz^2 -
u(i,j,k)*(2/hx^2 + 2/hy^2 + 2/hz^2) = c(x,y,z)
for i=1,nx-2 j=1,ny-2 k=1,nz-2 while substituting
u(i,j,k) boundary values when i=0 or nx-1, j=0 or ny-1,
k=0 or nz-1 using x=xmin+i*hx, y=ymin+j*hz, z=zmin+k*hz
Remember: hx, hy, hz are known constants, c(x,y,z) is computable.
The matrix has m = (nx-2)*(ny-2)*(nz-2) rows and columns for
solving for m equations in m unknowns. The right hand side vector
is the c(x,y,z) values with possibly boundary terms subtracted.
The unknowns are u(i,j,k) for i=1,nx-2, j=1,ny-2, k=1,nz-2 .
The first row of the matrix is for the unknown u(1,1,1) i=1, j=1, k=1
ub(xmin+0*hx,ymin+1*hy,zmin+1*hz)/hx^2 + u(2,1,1)/hx^2 +
ub(xmin+1*hx,ymin+0*hy,zmin+1*hz)/hy^2 + u(1,2,1)/hy^2 +
ub(xmin+1*hx,ymin+1*hy,zmin+0*hz)/hz^2 + u(1,1,2)/hz^2 +
u(1,1,1)*(2/hx^2 + 2/hy^2 + 2/hz^2) = c(xmin+1*hx,ymin+1*hy,zmin+1*hz)
Matrix cells in first row are zero except for:
M(1,1,1)(1,1,1) = 1/(2/hx^2 + 2/hy^2 + 2/hz^2)
M(1,1,1)(1,1,2) = 1/hz^2
M(1,1,1)(1,2,1) = 1/hy^2
M(1,1,1)(2,1,1) = 1/hx^2
V(1,1,1) = c(xmin+1*hx,ymin+1*hy,zmin+1*hz)
-ub(xmin+0*hx,ymin+1*hy,zmin+1*hz)/hx^2
-ub(xmin+1*hx,ymin+0*hy,zmin+1*hz)/hy^2
-ub(xmin+1*hx,ymin+1*hy,zmin+0*hz)/hz^2
Note that three of the terms are boundary terms and thus have a
computed value that is subtracted from the right hand side vector V.
The row of the matirx for the u(2,2,2) unknown has no boundary terms
u(1,2,2)/hx^2 + u(3,2,2)/hx^2 +
u(2,1,2)/hy^2 + u(2,3,2)/hy^2 +
u(2,2,1)/hz^2 + u(2,2,3)/hz^2 +
u(2,2,2)*(2/hx^2 + 2/hy^2 + 2/hz^2) = c(xmin+2*hx,ymin+2*hy,zmin+2*hz)
Matrix cells in row (2,2,2) are zero except for:
M(2,2,2)(1,2,2) = 1/hx^2
M(2,2,2)(3,2,2) = 1/hx^2
M(2,2,2)(2,1,2) = 1/hy^2
M(2,2,2)(2,3,2) = 1/hy^2
M(2,2,2)(2,2,1) = 1/hz^2
M(2,2,2)(2,2,3) = 1/hz^2
M(2,2,2)(2,2,2) = 1/(2/hx^2 + 2/hy^2 + 2/hz^2)
V(2,2,2) = c(xmin+2*hx,ymin+2*hy,zmin+2*hz)
The matrix is quite sparse and may be considered band diagonal.
There are only four possible values in the matrix. These are
computed only once at the start of building the matrix.
Solving the system of equations for our test problem gave
the error on the order of 10^-12 compared with
the iterative solution error of 10^-4 after 150 iterations.
The same code has been programmed in "C", Java, Fortran 95, Ada 95 and
Matlab 7 as shown below with file types .c, .java, .f90, .adb and .m:
"C" source code pde3b_eq.c
output of pde3b_eq.c
Fortran 95 source code pde3b_eq.f90
output of pde3b_eq.f90
Java source code pde3b_eq.java
output of pde3b_eq.java
Ada 95 source code pde3b_eq.adb
output of pde3b_eq.adb
Matlab 7 source code pde3b_eq.m
output of pde3b_eq.m
Please note that the solution of partial differential equations
gets significantly more complicated if the solution is not
continuous or is not continuously differentiable.
You can tailor the code presented above for many second order partial
differential equations in three independent variables.
1) You must provide a function u(x,y,z) that computes the solution
on the boundaries. It does not have to be a complete solution and
the "check" code should be removed.
2) You must provide a function c(x,y,z) that computes the
forcing function, the right hand side of the differential equation.
3a) If the differential equation is not d^2/dx^2 + d^2/dy^2 +d^2/dz^2
then use coefficients from nderiv.out to create a new function
u000(i,j,k) that computes a new value of the solution at i,j,k
from the neighboring coordinates for the next iteration.
3b) If the differential equation is not d^2/dx^2 + d^2/dy^2 +d^2/dz^2
then use coefficients from nderiv.out to create a new matrix
initialization function and then solve the simultaneous equations.
For further reading, more advanced, see:
Partial Differential Equations
<- previous index next ->