<- previous index next ->
Some partial differential equations are difficult to solve symbolically
and some are difficult to solve numerically.
This lecture covers the numerical solution of simple partial
differential equations where the boundary values are known.
See simple definitions if needed.
We are told there is an unknown function z(x,y) that satisfies
the partial differential equation:
(read 'd' as the partial derivative symbol, '^' as exponent)
d^2z(x,y)/dx^2 + d^2z(x,y)/dy^2 = c(x,y)
and we are given the function c(x,y)
and we are given the values of z(x,y) on the boundary
of the domain of x,y in xmin,ymin to xmax,ymax
To compute the numerical solution of z(x,y) at a set of points
we choose the number of points in the x and y directions, nx, ny.
Thus we have a grid of cells typically programmed as a
two dimensional array.
From the above we get a step size hx = (xmax-xmin)/(nx-1) and
hy = (ymax-ymin)/(ny-1). For our first solution we choose
hx = hy = h and just use a step size of h.
If we need z(x,y) at some point not on our solution grid,
we can use two 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^2z(x,y)/dx^2 + d^2z(x,y)/dy^2 =
1/h^2(z(x-h,y)-2z(x,y)+z(x+h,y)) +
1/h^2(z(x,y-h)-2z(x,y)+z(x,y+h)) +O(h^2)
Setting the above equal to c(x,y) and collecting terms
d^2z(x,y)/dx^2 + d^2z(x,y)/dy^2 = c(x,y) becomes
1/h^2(z(x-h,y)+z(x+h,y)+z(x,y-h)+z(x,y+h)-4z(x,y))=c(x,y)
now solve for a new z(x,y), based on neighboring grid points
z(x,y)= (z(z-h,y)+z(x+h,0)+z(x,y-h)+z(x,y+h)-c(x,y)*h^2)/4
The formula for z(x,y) is programmed based on indexes i and j
such that x=xmin+i*h, y=ymin+j*h for i=1,nx-2 j=1,ny-2
Note that the boundary values are set and unchanged for the four
sides i=0 for j=0,ny-1 i=nx-1 for j=0,ny-1
j=0 for i=0,nx-1 j=ny-1 for i=0,nx-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 z(x,y) 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, Ada 95 and
Matlab as shown below with file types .c, .java, .f90, .adb and .m:
"C" source code pde3.c
output of pde3.c
Fortran 95 source code pde3.f90
output of pde3.f90
Java source code pde3.java
output of pde3.java
Ada 95 source code pde3.adb
output of pde3.adb
Matlab 7 source code pde3.m
output of pde3.m
You may tailor any of the above code to solve second order
partial differential equations in two independent variables.
1) You need a function c(x,y) that is the forcing function, the right
hand side of the differential equation.
2) You need a function z(x,y) that provides at least the boundary
conditions. If your z(x,y) is not a test solution, remove the "check"
code.
3) If your differential equation is different, then use nderiv.out
to find the coefficients and create the function z00(i,j) that
computes the next value of a solution point from the surrounding
previous solution points.
The following examples will be covered for a more general case
in the next lecture. The same problem is solved as given above,
using a method of solving simultaneous linear equations, rather
than an iterative method. The files corresponding to the above
examples are with minimum modifications are:
"C" source code pde3_eq.c
output of pde3_eq.c
Matlab 7 source code pde3_eq.m
output of pde3_eq.m
The instructor understands that some students have a strong
prejudice in favor of, or against, some programming languages.
After about 50 years of programming in about 50 programming
languages, the instructor finds that the difference between
programming languages is mostly syntactic sugar. Yet, since
students may be able to read some programming languages
easier than others, these examples are presented in "C",
Fortran 90, Java, Ada and Matlab. The intent was to do a quick
translation, keeping most of the source code the same,
for the different languages. Style was not a consideration.
Some rearranging of the order was used when convenient.
In Java and Ada the quick and dirty formatting was used
and thus it is suggested to look at the formatted output of
"C" or Fortran. The numerical results are almost exactly
the same although the presentation differs considerably.
The Matlab code is a combination of "C" and Fortran given
as a single file with internal functions. The subscript base
is shifted from 0 to 1 as in old Fortran.
<- previous index next ->