# Help with NDSolveFiniteDifferenceDerivative

Posted 9 years ago
8483 Views
|
6 Replies
|
8 Total Likes
|
 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 NDSolveFiniteDifferenceDerivative, 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.
6 Replies
Sort By:
Posted 9 years ago
 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 9 years ago
 Thanks, Daniel. Is there a simple way to map your output to an fdd matrix?
Posted 9 years ago
 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 9 years ago
 Thanks, Daniel -- that's perfect and solves my problem.
Posted 9 years ago
 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 = NDSolveFiniteDifferenceDerivative[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 9 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, ..., 25Because 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?