<- previous    index    next ->

Lecture 31, Creating PDE Test Cases

The message here:
Thoroughly test any software you write or are about to use.

Use tools to generate test when possible.

Generate a test case for a second order PDE with two independent variables.
Make the PDE general by having it be  elliptic in some regions,
hyperbolic in some regions and passing through parabolic.

Note: This is for the set of solutions that are continuous and
continuously differentiable. There are many special solutions and
corresponding special solvers that should be used for special cases.
This lecture is about a general case and a general solver.

The solution will be u(x,y).
The notation for derivatives will be
  uxx(x,y) second derivative of u(x,y) with respect to x
  uxy(x,y) derivative of u(x,y) with respect to x and with respect to y
  uyy(x,y) second derivative of u(x,y) with respect to 
  ux(x,y) derivative of u(x,y) with respect to x
  uy(x,y) derivative of u(x,y) with respect to y

We will create functions a1(x,y), b1(x,y), c1(x,y), d1(x,y), e1(x,y) and
f1(x,y).

The PDE we intend to solve, compactly, is:

 a1 uxx + b1 uxy + c1 uyy + d1 ux + e1 uy + f1 u = c

The PDE written for computer processing is:

 a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
  d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y)*u(x,y) = c(x,y)

Yes, we could use algebra to reorganize the PDE, yet we will not.

We choose the functions a1, b1 and c1 such that b1^2 -a1 * c1 is
elliptic, hyperbolic and parabolic in various regions.

 a1(x,y) = exp(x/2)*exp(y)/2

 b1(x,y) = 0.7/(x*x*y*y + 0.5)

 c1(x,y) = (4 - exp(x) - exp(y/2))*2

then we plot  b1^2 - a1 * c1 over the x-y plane.

abc.m MatLab to plot a1, b1, c1 and b1^2-a1*c1






Below the yellow plane is elliptic, above the yellow plane is
hyperbolic and on the yellow plane is parabolic. This has little
to do with the numerical solution of our PDE. For some classic
PDE's only in one region, there are specialized solutions.


Equations of mixed type:
If a PDE has coefficients which are not constant, as shown above,
it is possible that it will not belong to any of the three categories.
A simple but important example is the Euler-Tricomi equation
    uxx = xuyy
which is called elliptic-hyperbolic because it is elliptic in the region
x < 0, hyperbolic in the region x > 0, and degenerate parabolic
on the line x = 0.


The other functions for this test case were chosen as:

 d1(x,y) = x^2 y

 e1(x,y) = x y^2

 f1(x,y) = 3x + 2y


Now we have to pick our solution function u(x,y).
Since we will be taking second derivatives, we need at least a
third order solution in order to be interesting.

Out of the air, I chose:

  u(x,y) = x^3 + 2y^3 + 3x^2 y + 4xy^2 + 5xy + 6x + 7y + 8

We do not want the solution symmetric in x and y.
We want a smattering of terms for testing.
I use 1, 2, 3, 4, ... to be sure each term tested.
Now we have to compute the forcing function   c(x,y)
Doing the computation by hand is too error prone.
I use Maple.

Maple computation of c(x,y)  Look at this!

Maple work sheet
Maple output

c(x,y) =  0.5*exp(x/2.0)*exp(y)*(6.0*x+6.0*y)                  +
          0.7*(6.0*x + 8.0*y + 5.0)/(x*x*y*y+0.5)              +
          (8.0 - 2.0*exp(x) - 2.0*exp(y/2.0))*(12.0*y + 8.0*x) +
          (x*x+y)*(3.0*x*x + 6.0*x*y + 4.0*y*y + 5.0*y + 6.0)  +
          x*y*y*(6.0*y*y + 3.0*x*x + 8.0*x*y + 5.0*x +7.0)     +
          (3.0*x + 2.0*y)*(x^3 + 2y^3 + 3x^2 y + 4xy^2 + 5xy + 6x + 7y + 8)


We have to code the function  c(x,y) for our solver.
We have to code the function  u(x,y) for our solver to compute
boundary values and we will use the function to check our solver.

The solution will be u(x,y) and the solver must be able to evaluate
this function on the boundary. The region will be
  xmin <= x <= xmax
  ymin <= y <= ymax
We will solve for nx points in the x dimension, ny points in y dimension with 
  hx=(xmax-xmin)/(nx-1)  step size in x
  hy=(ymax-ymin)/(ny-1)  step size in y

We will compute the approximate solution at u(i,j) for the point
u(xmin+i*hx,ymin+j*hy). With subscripts running i=0..nx-1 and j=0..ny-1.
Known boundary values are at i=0, i=nx-1, j=0, and  j=ny-1.
These subscripts are for "C" and Java, add one for Ada, Fortran and
MatLab subscripts.


The actual numeric values can be set at the last minute, yet we will
use  xmin = ymin = -1.0    xmax = ymax = 1.0   nx = 7   ny = 7

Now, we have to generate the matrix that represents the system of
linear simultaneous equations for the unknown values of u(x,y)
at xmin+i*hx  for i=1..nx-2    ymin+j*hy  for j=1..ny-2

I am using solution points  1, 2, ... , nx-2 by
1, 2, ... , ny-2 stored in a matrix starting at ut[0][0] for coding in "C".

The solution points will be the same as used in pde3:



The PDE will be made discrete by using unknown values in difference
equations as given in nderiv.out or
computed by deriv.c, deriv.adb, deriv.java, deriv.py, deriv.scala, etc.

For this test case, I choose to use five point derivatives for
both first and second order. Note that the uxy term has the
first derivative with respect to x at five y values then
the first derivative of these values with respect to y.


Taking a term at a time from the PDE and writing the discrete version:

 a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
  d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y)*u(x,y) = c(x,y)

At (x,y)  the discrete approximation of  a1(x,y)*uxx(x,y) =

 a1(x,y) * (1/(12*hx*hx)) *
 (-1*u(x-2hx,y) +16*u(x-hx,y) -30*u(x,y) +16*u(x+hx,y) -1*u(x+2hx,y))

yet, we want to solve for u(i,j) in a matrix,
thus using x=xmin+i*hx and y=ymin+j*hy rewrite the above approximation as

 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hx)) *
 (-1*u(i-2,j) +16*u(i-1,j) -30*u(i,j) +16*u(i+1,j) -1*u(i+2,j))


At (x,y)  c1(x,y)*uyy(x,y) = the discrete approximation

 c1(x,y) * (1/(12*hy*hy)) *
 (-1*u(x,y-2hy) +16*u(x,y-hy) -30*u(x,y) +16*u(x,y+hy) -1*u(x,y+2hy))

again, we want to solve for u(i,j) in a matrix,
thus using x=xmin+i*hx and y=ymin+j*hy rewrite the above approximation as

 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *
 (-1*u(i,j-2) +16*u(i,j-1) -30*u(i,j) +16*u(i,j+1) -1*u(i,j+2))


What happened to b1(x,y)*uxy(x,y) ?
Well that one is more difficult, thus we now use the combination of
the first derivative with respect to x, then the first derivative
with respect to y, 25 terms.
Writing all the terms using i and j

b1(xmin+i*hx,ymin+j*hy) * (1/(144*hx*hy)) *
( 1*u(i-2,j-2)  -8*u(i-1,j-2) +0*u(i,j-2)  +8*u(i+1,j-2) -1*u(i+2,j-2)
 -8*u(i-2,j-1) +64*u(i-1,j-1) +0*u(i,j-1) -64*u(i+1,j-1) +8*u(i+2,j-1)
 +0*u(i-2,j  )  +0*u(i-1,j  ) +0*u(i,j  )  +0*u(i+1,j  ) +0*u(i+2,j  )
 +8*u(i-2,j+1) -64*u(i-1,j+1) +0*u(i,j+1) +64*u(i+1,j+1) -8*u(i+2,j+1)
 -1*u(i-2,j+2)  +8*u(i-1,j+2) +0*u(i,j+2)  -8*u(i+1,j+2) +1*u(i+2,j+2))

Just writing the remaining terms using i and j
using the discrete coefficients.

  d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *
  (1*u(i-2,j) -8*u(i-1,j) +0*u(i,j) +8*u(i+1,j) -1*u(i+2,j))

  e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *
  (1*u(i,j-2) -8*u(i,j-1) +0*u(i,j) +8*u(i,j+1) -1*u(i,j+2))


The u(x,y) term is just

  f1(xmin+i*hx,ymin+j*hy)*u(i,j)

The right hand side of every approximation above is just

  c(xmin+i*hx,ymin+j*hy)


If any of the u(i+a,j+b) are boundary elements, we plug in the
numeric value of the boundary element, multiply by the
coefficient, and subtract the product from the right hand side.
 
Now we have many equations for u(i,j) using values for i and j,
and we must compute the coefficients for the set of equations
for  i=2..nx-3  j=2..ny-3  the central case.

Remember, u(i=0,j), u(i=nx-1,j), u(i,j=0) and u(i,j=ny-1) are
known boundary values.

What about u(1,j), u(nx-2,j) for j=1..ny-2 and
           u(i,1), u(i,ny-2) for i=2..nx-3  (do not use u(1,1) twice !!!)
                                            (do not use u(nx-2,ny-2) twice)

Fortunately, the coefficients for the discrete derivatives can be
computed by deriv.c, deriv.adb, etc.

At (1,j) the discrete approximation of 
         a1(xmin+1*hx,ymin*j*hy)*uxx(xmin+1*hx,ymin+j*hy) =

 a1(xmin+1*hx,ymin+j*hy) * (1/(12*hx*hx)) *
 (-3*u(i-1,j) -10*u(i,j) +18*u(i+1,j) -6*u(i+2,j) +1*u(i+3,j))

 We could not use i-2 because it is outside our region, thus
 note "deriv" is called with 'point' set to 1.


At (nx-2,j)  the discrete approximation of
             a1(xmin+(nx-2)*hx,ymin*j*hy)*uxx(xmin+(nx-2)*hx,ymin+j*hy) =

 a1(xmin+(nx-2)*hx,ymin+j*hy) * (1/(12*hx*hx)) *
 (-1*u(i-3,j) +6*u(i-2,j) -18*u(i-1,j) +10*u(i,j) +3*u(i+1,j))

 We could not use i+2 or (nx-2)+2 because it is outside our region, thus
 note "deriv" is called with 'point' set to 3.

There is nothing special about uxx(i,1) or uxx(nx-2,1) use the general case.



At (i,1)  the discrete approximation of  c1(x,y)*uyy(x,y) =

 c1(xmin+i*hx,ymin+1*hy) * (1/(12*hy*hy)) *
 (-3*u(i,j-1) -10*u(i,j) +18*u(i,j+1) -6*u(i,j+2) +1*u(i,j+3))

 We could not use j-2 because it is outside our region, thus
 note "deriv" is called with 'point' set to 1.


At (i,nx-2)  the discrete approximation of  c1(x,y)*uyy(x,y) =

 c1(xmin+i*hx,ymin+(nx-2)*hy) * (1/(12*hy*hy)) *
 (-1*u(i,j-3) +6*u(i,j-2) -18*u(i,j-1) +10*u(i,j) +3*u(i,j+1))

 We could not use j+2 or (ny-2)+2 because it is outside our region, thus
 note "deriv" is called with 'point' set to 3.

There is nothing special about uyy(1,j) or uxx(nx-2,j) use the general case.

Oh, and yes, ux(x,y) uy(x,y) uxy(x,y) also have to be shifted for the
"just inside the boundary case." Think of how much easier it is using five
points rather than sever points. With seven points the two rows and
columns inside the boundary are special cases.


Now we have enough equations to exactly compute the approximate solution.
We build a system of linear equations of the form:

  | ut00  ut01  ut02  ut03 |   | u(1,1) |   | k0 |
  | ut10  ut11  ut12  ut13 | * | u(1,2) | = | k1 |
  | ut20  ut21  ut22  ut23 |   | u(2,1) |   | k2 |
  | ut30  ut31  ut32  ut33 |   | u(2,2) |   | k3 |

We know the values at the boundary u(0,0), u(0,1), u(0,2), u(0,3)
                                   u(1,0),                 u(1,3)
                                   u(2,0),                 u(2,3)
                                   u(3,0), u(3,1), u(3,2), u(3,3)
 
For this specific system of equations, nx=4, ny=4 and there are four
internal, non boundary, values to be found. The number of equations
will always be (nx-2)*(ny-2).

The value found for u(1,1) from solving the system of linear equations
is the desired solution at u(xmin+1*hx,ymin+1*hy). u(2,2) is the value of
u(xmin+2*hx,ymin+2*hy). Additional values, not on hx or hy steps,
of u(x,y) may be found by two dimensional interpolation.

The matrix "ut", used in the source code pde_abc_eq.c
is cleared to zero. Then each equation is used and the
coefficient of u(i,j) is added the appropriate ut entry.

There will be (nx-2)*(ny-2) equations that must be used.
The i and j for these equations are  for i=1..nx-2  for j=1..ny-2.
The PDE approximation above must be used for these i,j pairs.
Note that u(i,j), u(i-2,j) etc are variables and it is the
coefficients of these variables that get added into the "ut" matrix.
The first equation will add to entries in the first row of the matrix.
The second equation will add to entries in the second row of the matrix.
The boundary values are not stored in the ut matrix and using algebra,
a boundary value coefficient is multiplied by the boundary value and
subtracted from the constant term. This can cause many special cases
in the solver.
In general many entries in a row will not be changed from zero.
The final matrix will be in the class of a band matrix.

Now the gruesome work of evaluating the u(i,j) coefficients.
We cheat here. Now assume nx>4 and ny>5 and thus
i=2, j=3 is a general center case. From above, we have


 a1(x,y)*uxx(x,y) + b1(x,y)*uxy(x,y) + c1(x,y)*uyy(x,y) +
  d1(x,y)*ux(x,y) + e1(x,y)*uy(x,y)  + f1(x,y)*u(x,y) = c(x,y)

converted to a discrete approximation for the central elements is

 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hx)) *
  (-1*u(i-2,j) +16*u(i-1,j) -30*u(i,j) +16*u(i+1,j) -1*u(i+2,j))          +
 b1(xmin+i*hx,y+j*hy) * (1/(144*hx*hy)) *
  ( 1*u(i-2,j-2)  -8*u(i-1,j-2) +0*u(i,j-2)  +8*u(i+1,j-2) -1*u(i+2,j-2)
   -8*u(i-2,j-1) +64*u(i-1,j-1) +0*u(i,j-1) -64*u(i+1,j-1) +8*u(i+2,j-1)
   +0*u(i-2,j  )  +0*u(i-1,j  ) +0*u(i,j  )  +0*u(i+1,j  ) +0*u(i+2,j  )
   +8*u(i-2,j+1) -64*u(i-1,j+1) +0*u(i,j+1) +64*u(i+1,j+1) -8*u(i+2,j+1)
   -1*u(i-2,j+2)  +8*u(i-1,j+2) +0*u(i,j+2)  -8*u(i+1,j+2) +1*u(i+2,j+2)) +
 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *
  (-1*u(i,j-2) +16*u(i,j-1) -30*u(i,j) +16*u(i,j+1) -1*u(i,j+2))          +
 d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *
  (1*u(i-2,j) -8*u(i-1,j) +0*u(i,j) +8*u(i+1,j) -1*u(i+2,j))              +
 e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *
  (1*u(i,j-2) -8*u(i,j-1) +0*u(i,j) +8*u(i,j+1) -1*u(i,j+2))              +
 f1(xmin+i*hx,ymin+j*hy)*u(i,j)                                           =
 c(xmin+i*hx,ymin+j*hy)

Collecting terms for u(i,j) we compute the coefficient  ct
 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hx)) *(-30) +
 b1(xmin+i*hx,ymin+j*hy) * (1/(144*hx*hy)) *(0)     +
 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *(-30) +
 d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *(0)       +
 e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *(0)

and the right hand side is  c(xmin+i*hx,ymin+j*hy)

A bit of bizarre subscripting, assuming we are at i=2, j=3 this
is equation number ii=(i-1)+(nx-2)*(j-1), right?
Check i=1,j=1    is equation 0    if i or j equal zero, we have a boundary
      i=2,j=1    is equation 1
      ...
      i=nx-2,j=1 is equation nx-3
      i=1,j=2    is equation nx-2
      i=2,j=2    is equation nx-1
      ...
      i=2,j=3    is equation (2-1)+(nx-2)*(3-1) the row in the ut matrix.

cs, the subscript for the right hand side is  cs=(nx-2)*(ny-2) stored in ut.

Thus we add ct to ut(ii,ii). The second subscript is where u(i,j) is.
Well, to tell the truth, we add each component of the sum in a loop,
directly into ut(ii,ii). If we had  ct  we would just store it in ut(ii,ii).

ut(ii,cs) is set to  c(xmin+i*hx,ymin+j*hy) - f1(xmin+i*hx,ymin+j*hy)


We must, of course, compute all the coefficients, e.g. u(i-1,j), ctm1

 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hx)) *(16) +
 b1(xmin+i*hx,ymin+j*hy) * (1/(144*hx*hy)) *(0)    +
 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *(0)  +  does not exists
 d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *(-8)     +
 e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *(0)         does not exists

ij=(((i-1)-1)+(nx-2)*(j-1) where u(i-1,j) is in the solution vector.
ut(ii,ij) = ctm1

We must, of course, compute all the coefficients, e.g. u(i-2,j), ctmm1

 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hy)) *(-1) +
 b1(xmin+i*hx,ymin+j*hy) * (1/(144*hx*hy)) *(0)    +
 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *(0)  +  does not exists
 d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *(1)     +
 e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *(0)         does not exists

ij=(((i-2)-1)+(nx-2)*(j-1) where u(i-2,j) is in the solution vector.
ut(ii,ij) = ctmm1   Wrong if i is 2!  u(0,j) is a boundary value, thus
ut(ii,cs) = ut(ii,cs) - ctmm1 * u(0,j) 

And, one more coefficients, e.g. u(i,j+1), ctp1

 a1(xmin+i*hx,ymin+j*hy) * (1/(12*hx*hx)) *(0)  +  does not exists
 b1(xmin+i*hx,ymin+j*hy) * (1/(144*hx*hy)) *(0)    +
 c1(xmin+i*hx,ymin+j*hy) * (1/(12*hy*hy)) *(16) +
 d1(xmin+i*hx,ymin+j*hy) * (1/(12*hx) *(0 )     +  does not exists
 e1(xmin+i*hx,ymin+j*hy) * (1/(12*hy) *(8)

"C" and Java subscripting
ij=((i-1)+(nx-2)*((j+1)-1) when u(i,j+1) is in the solution vector.
ut(ii,ij) = ctp1
Ada, Fortran and MatLab subscripting
ij=(i-1)+(nx-2)*((j+1)-2) when u(i,j+1) is in the solution vector.
ut(ii,ij) = ctp1
where uc(ii) is the right hand side.

Yes, the code is nested four levels deep in iterations. 

Now do u(1,j), u(nx-2,j) for j=1..ny-2 and
       u(i,1), u(i,ny-2) for i=2..nx-3  (do not use u(1,1) twice !!!)
                                        (do not use u(nx-2,ny-2) twice)


Once the "ut" matrix is initialized, solve the simultaneous equations
and print the answers. Now, wasn't that easy.

A crude layout for  nx=9, ny=9 is shown below with subscripts
The  i,j  notation is for a known boundary value,
the number is the ii subscript of the unknown value inside the I====I
                                                               I ut I
                                                               I====I.

                               j
       0     1     2     3     4     5     6     7     8   ny=9
    +-----------------------------------------------------+
  0 | 0,0 | 0,1 | 0,2 | 0,3 | 0,4 | 0,5 | 0,6 | 0,7 | 0,8 |
    |     I=========================================I     |
  1 | 1,0 I  0  |   7 |  14 |  21 |  28 |  35 |  42 I 1,8 |
    |     I-----------------------------------------I     |
  2 | 2,0 I  1  |   8 |  15 |  22 |  29 |  36 |  43 I 2,8 |
    |     I-----------------------------------------I     |
  3 | 3,0 I  2  |   9 |  16 |  23 |  30 |  37 |  44 I 3,8 |
    |     I-----------------------------------------I     |
i 4 | 4,0 I  3  |  10 |  17 |  24 |  31 |  38 |  45 I 4,8 |
    |     I-----------------------------------------I     |
  5 | 5,0 I  4  |  11 |  18 |  25 |  32 |  39 |  46 I 5,8 |
    |     I-----------------------------------------I     |
  6 | 6,0 I  5  |  12 |  19 |  26 |  33 |  40 |  47 I 6,8 |
    |     I-----------------------------------------I     |
  7 | 7,0 I  6  |  13 |  20 |  27 |  34 |  41 |  48 I 7,8 |
    |     I=========================================I     |
  8 | 8,0 | 8,1 | 8,2 | 8,3 | 8,4 | 8,5 | 8,6 | 8,7 | 8,8 |
    +-----------------------------------------------------+
 nx=9
           ut(ii) where ii = (i-1)+(nx-2)*(j-1)
These are "C" and Java subscripts.

                               j
       1     2     3     4     5     6     7     8     9   ny=9
    +-----------------------------------------------------+
  1 | 1,1 | 1,2 | 1,3 | 1,4 | 1,5 | 1,6 | 1,7 | 1,8 | 1,9 |
    |     I=========================================I     |
  2 | 2,1 I  1  |   8 |  15 |  22 |  29 |  36 |  43 I 2,9 |
    |     I-----------------------------------------I     |
  3 | 3,1 I  2  |   9 |  16 |  23 |  30 |  37 |  44 I 3,9 |
    |     I-----------------------------------------I     |
  4 | 4,1 I  3  |  10 |  17 |  24 |  31 |  38 |  45 I 4,9 |
    |     I-----------------------------------------I     |
i 5 | 5,1 I  4  |  11 |  18 |  25 |  32 |  39 |  46 I 5,9 |
    |     I-----------------------------------------I     |
  6 | 6,1 I  5  |  12 |  19 |  26 |  33 |  40 |  47 I 6,9 |
    |     I-----------------------------------------I     |
  7 | 7,1 I  6  |  13 |  20 |  27 |  34 |  41 |  48 I 7,9 |
    |     I-----------------------------------------I     |
  8 | 8,1 I  7  |  14 |  21 |  28 |  35 |  42 |  49 I 8,9 |
    |     I=========================================I     |
  9 | 9,1 | 9,2 | 9,3 | 9,4 | 9,5 | 9,6 | 9,7 | 9,8 | 9,9 |
    +-----------------------------------------------------+
 nx=9
           ut(ii) where ii = (i-1)+(nx-2)*(j-2)
These are Ada, Fortran and Matlab subscripts.


A differential equation that has regions hyperbolic, elliptic and parabolic
defined in the various languages "abc" file.

The basic functions coded in "C" for this test case are:

abc.txt  instructions for user
abc.h  sample test case described above
abc.c  implementation of test case
deriv.h  discretization header file
deriv.c  discretization
pde_abc_eq.c  solver program
pde_abc_eq_c.out  sample output
abc.c  implementation of test case

The basic functions coded in Fortran 90 for this test case are:

abc.f90  implementation of test case
deriv.f90  discretization
simeq.f90  solve simultaneous equations
pde_abc_eq.f90  solver program
pde_abc_eq_f90.out  sample output

The basic functions coded in Java for this test case are:

abc.java  sample test case described above
nderiv.java  discretization
simeq.java  solve simultaneous equations
pde_abc_eq.java  solver program
pde_abc_eq_java.out  sample output

The basic functions coded in Ada 95 for this test case are:

abc.ads  sample test case described above
abc.adb  implementation of test case
deriv.adb  discretization
rderiv.adb  discretization
simeq.adb  solve simultaneous equations
pde_abc_eq.adb  solver program
pde_abc_eq_ada.out  sample output



Sparse matrix storage is needed for the system of linear equations
when nx or ny is significantly greater than the discretization
order, nd.

The same PDE's as above, using the definition in "abc" are
easily modified to use sparse matrix storage and a direct
sparse matrix solution to the linear equations.

The basic functions coded in "C" for this test case are:

abc.txt  instructions for user
sparse.h  sparse matrix header file
sparse.c  sparse matrix for PDE's
deriv.h  discretization header file
deriv.c  discretization
abc.h  sample test case described above
abc.c  implementation of test case
sparse_abc.c  sparse solver program
sparse_abc_c.out  sample output

The basic functions coded in Java for this test case are:

abc.java  sample test case described above
sparse.java  sparse matrix for PDE's
nderiv.java  discretization
sparse_abc.java  solver program
sparse_abc_java.out  sample output

The basic functions coded in Ada 95 for this test case are:

abc.ads  sample test case described above
abc.adb  implementation of test case
sparse.ads  sparse matrix specification
sparse.adb  sparse matrix for PDE's
deriv.adb  discretization
sparse_abc.adb  solver program
sparse_abc_ada.out  sample output




Some challenging fourth order partial differential equations
including biharmonic equations 
are covered in lecture 28d

Solving nonlinear PDE is covered in Lecture 31b.
The specific case for the nonlinear Navier Stokes Equations is
covered in Lecture 28b

For nonlinear PDE to be solved by solving a system of equations,
the simultaneous equations are non linear.
To solve a system of nonlinear equations reasonably efficiently,
use a Jacobian matrix and an iterative solution.

An example that can handle up to third order nonlinear systems of
equations is shown by using a Newton iteration based on the inverse
of the Jacobian:

  x_next = x_prev - J^-1 * ( A * x_prev - Y)

simeq_newton3.java basic solver
test_simeq_newton3.java test program
test_simeq_newton3_java.out test_output
on three test cases (minor heuristic used to build J)



Plotting Solutions

Plotting solutions can be easy with Linux and gnuplot. To demonstrate, I have a program that generates a 2D solution that needs a 3D plot. make_gnuplotdata.c The data file, for this case, is gnuplot.dat The input file to gnuplot is my2.plot The shell script to make the plot is test2_gnuplot.sh The plot is run by typing sh test2_gnuplot.sh or from inside a program with test2_gnuplot.c or from command line gnuplot -persist my2.plot The result is the plot

A special case is tests for Laplace Equation

ΔU = 0 or ∇2U = 0 or Uxx(x,y) + Uyy(x,y) = 0 or ∂2U/∂x2 + ∂2U/∂y2 = 0 Many known solutions and thus Dirichlet boundary conditions for test cases are known: k>0 is any constant U(x,y)=k U(x,y)=k*x U(x,y)=k*y U(x,y)=k*x*y U(x,y)=k*x^2 - k*y^2 U(x,y)=sin(k*x)*sinh(k*y) U(x,y)=sin(k*x)*cosh(k*y) U(x,y)=cos(k*x)*sinh(k*y) U(x,y)=cos(k*x)*cosh(k*y) U(x,y)=sin(k*x)*exp( k*y) U(x,y)=sin(k*x)*exp(-k*y) U(x,y)=cos(k*x)*exp( k*y) U(x,y)=cos(k*x)*exp(-k*y) r=sqrt(x^2+y^2) th=atan2(y,x) -Pi..Pi principal value is U(x,y)=r^k * sin(k*th) discontinuous along negative X axis using atan2 in 0..2Pi wrong! causes discontinuous derivative along positive X axis
    <- previous    index    next ->

Other links

Go to top