Message Boards Message Boards

2
|
10177 Views
|
6 Replies
|
8 Total Likes
View groups...
Share
Share this post:

Help with NDSolve`FiniteDifferenceDerivative

Alan Lewis
Posted 11 years ago

Consider the operator $f_{xx} + f_{yy}$ in the unit square with f=0 on the boundaries. I would like to construct the differentiation matrix associated to that on a uniform tensor product grid in the interior of the square with the usual O(h^2) centered finite difference formula (5 point stencil).

If I apply NDSolve`FiniteDifferenceDerivative, say with difference order 2, as discussed in the NDSolve tutorial, I will get a matrix with O(h^2) interior pointing derivatives near the boundaries. I don't want that -- I simply want the usual centered differences. For example, if my xgrid runs x = (h, 2h, 3h, ...), the matrix should represent the difference formula (showing the x-coordinates only)

$f_{xx}(h) = \frac{f(2h) - 2 f(h) + f(0)}{h^2} = \frac{f(2h) - 2 f(h)}{h^2}$ at that first grid point.

Is there a way to use the NDSolve utilities to get this? After playing around, I couldn't see how, and ended up constructing my differentiation matrix "by hand", in a For loop running through all the grid points.
POSTED BY: Alan Lewis
6 Replies
Sort By:
Alan Lewis
Alan Lewis
Posted 11 years ago

Thanks, Daniel -- that's perfect and solves my problem.
POSTED BY: Alan Lewis

Not sure if I understand but maybe you want something like this?

(* Repeat setup sand name the resulting 5x5 *)
kernel = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}};
vals = Array[f, {5, 5}];
laplacian = ListConvolve[kernel, vals, {2, -2}, 0];

Now we regard this as 25 expressions in 25 variables and extract the coefficient arrays. The first is the (missing) constants, which would be the right hand side in an inhomogeneous equation of the form mat.x==b. We ignore it here (it's all zeros).

lapOperator = CoefficientArrays[Flatten@laplacian, Flatten@vals][[2]]

(* Out[571]= SparseArray[< 105 >, {25, 25}] *)

The sparse array form is probably preferable for actual use, but if you want to see what it looks like here it is in full matrix form.

Normal[lapOperator]

(* Out[572]= {{-4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0}, {1, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 1, -4, 1, 0, 0, 0, 1,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, -4, 
  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {1, 0, 
  0, 0, 0, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0}, {0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0}, {0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 
  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 1, 0, 
  0, 0, 0, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0,
   0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0,
   0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 
  0, 1, 0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 
  1, -4, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   1, 0, 0, 0, 0, -4, 1, 0, 0, 0, 1, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 0, 0}, {0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0, 0, 1, 0, 
  0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0,
   0, 0, 1, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0,
   1, -4, 0, 0, 0, 0, 1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 
  0, 1, 0, 0, 0, 0, -4, 1, 0, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
   0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0, 0}, {0, 0, 0, 0, 0, 0, 0, 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 1, 0}, {0, 0, 0, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, -4, 
  1}, {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 
  0, 0, 1, -4}} *)
POSTED BY: Daniel Lichtblau
Alan Lewis
Alan Lewis
Posted 11 years ago

Thanks, Daniel. Is there a simple way to map your output to an fdd matrix?
POSTED BY: Alan Lewis
Alan Lewis
Alan Lewis
Posted 11 years ago

Thanks, Rob.

I tried your scheme but it doesn't seem to do what I need. Let's take {nx,ny} = {5,5} for example with lattice cords x(i) = i h, y(j) = j h, (i,j)=1,...n for the interior of the unit square; h=1/(n+1)

The fdd has size n^2 x n^2 = 25 x 25. Let's number the fdd rows with the index IBIG = 1, 2, ..., 25

Because of the 5 point stencil, if a point in the square with corresponding ffd coord IBIG has 4 neighbors also in the interior, then the corresponding fdd row (with row number IBIG), should have 5 non-zero entries. (itself + neighbors)

However, if a point in the square has fewer than 4 nearest neigbors (because some of its neighbors have x=0 or y=0 or x=1 or y=1), then the corresponding fdd row should have fewer than 5 entries. For example, the first row of the fdd, which represents the stencil associated to the interior point in the lower left interior corner in (x,y) plane, should have only 3 non-zero entries.

Your fdd matrix, on the other hand, has 5 non-zero entries for every row ...

I appreciate the help. Can you see a fix?
POSTED BY: Alan Lewis

Rob Knapp showed how this might be done with NDSolve utility functions. The ListConvolve approach below might also be useful in cases where these utilities might be hard to fathom. The kernel matrix gives the stencil and coefficients for that difference operator. I'll illustrate on a 5x5 matrix of unknowns. I believe I got the overhangs correct for the stated purpose, but if not, adjust as needed.

kernel = {{0, 1, 0}, {1, -4, 1}, {0, 1, 0}};
vals = Array[f, {5, 5}];
ListConvolve[kernel, vals, {2, -2}, 0]

(* Out[436]= {{-4 f[1, 1] + f[1, 2] + f[2, 1], 
  f[1, 1] - 4 f[1, 2] + f[1, 3] + f[2, 2], 
  f[1, 2] - 4 f[1, 3] + f[1, 4] + f[2, 3], 
  f[1, 3] - 4 f[1, 4] + f[1, 5] + f[2, 4], 
  f[1, 4] - 4 f[1, 5] + f[2, 5]}, {f[1, 1] - 4 f[2, 1] + f[2, 2] + 
   f[3, 1], f[1, 2] + f[2, 1] - 4 f[2, 2] + f[2, 3] + f[3, 2], 
  f[1, 3] + f[2, 2] - 4 f[2, 3] + f[2, 4] + f[3, 3], 
  f[1, 4] + f[2, 3] - 4 f[2, 4] + f[2, 5] + f[3, 4], 
  f[1, 5] + f[2, 4] - 4 f[2, 5] + f[3, 5]}, {f[2, 1] - 4 f[3, 1] + 
   f[3, 2] + f[4, 1], 
  f[2, 2] + f[3, 1] - 4 f[3, 2] + f[3, 3] + f[4, 2], 
  f[2, 3] + f[3, 2] - 4 f[3, 3] + f[3, 4] + f[4, 3], 
  f[2, 4] + f[3, 3] - 4 f[3, 4] + f[3, 5] + f[4, 4], 
  f[2, 5] + f[3, 4] - 4 f[3, 5] + f[4, 5]}, {f[3, 1] - 4 f[4, 1] + 
   f[4, 2] + f[5, 1], 
  f[3, 2] + f[4, 1] - 4 f[4, 2] + f[4, 3] + f[5, 2], 
  f[3, 3] + f[4, 2] - 4 f[4, 3] + f[4, 4] + f[5, 3], 
  f[3, 4] + f[4, 3] - 4 f[4, 4] + f[4, 5] + f[5, 4], 
  f[3, 5] + f[4, 4] - 4 f[4, 5] + f[5, 5]}, {f[4, 1] - 4 f[5, 1] + 
   f[5, 2], f[4, 2] + f[5, 1] - 4 f[5, 2] + f[5, 3], 
  f[4, 3] + f[5, 2] - 4 f[5, 3] + f[5, 4], 
  f[4, 4] + f[5, 3] - 4 f[5, 4] + f[5, 5], 
  f[4, 5] + f[5, 4] - 4 f[5, 5]}} *)
POSTED BY: Daniel Lichtblau

You can use (for whatever dimension you like) where for 2D, h = {hx,hy}, n = {nx,ny}

zeroBoundaryLaplaceMatrix[h_List, n_List] /;  Length[h] == Length[n] := 
Module[{d = Length[h], fdd, coords},
    coords = MapThread[(#1 Range[0, #2]) &, {h, n}];
    Total[Table[
        fdd = NDSolve`FiniteDifferenceDerivative[2 UnitVector[d, i], coords, "DifferenceOrder" -> 2, "PeriodicInterpolation" -> True];
        fdd["DifferentiationMatrix"], {i, d}]]]

The zero boundary points in each dimension are at 0 and h*(n + 1)
POSTED BY: Rob Knapp
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract