Any partial differential equation can be converted from the continuous analytic form to a discrete form. I call the process discretization. The two dimensional example below is a fourth order PDE in x and y over a non uniform grid in both x and y. The unknown solution is U(x,y). The x points are given as xg[0]..xg[nx-1] and the y points are given as yg[0]..yg[ny-1]. Enough Dirichlet boundary conditions, k of them, are given at vertices (xg[ik],yg[jk]) to make the solution unique. The boundary vertices may be identified in many ways, for this example, an array ig[][] of size nx by ny has an indication of a bound or free vertex. The number of free vertices is the Degree of Freedom, DOF, of the specific problem. The, non trivial, example equation to be solved numerically is: Uxxxx(x,y) + 2 Uxxyy(x,y) + Uyyyy(x,y) + 2 Uxx(x,y) + 2 Uyy(x,y) + U(x,y) = f(x,y) Uxxxx(x,y) is the fourth derivative of U(x,y) with respect to x. The minimum discrete form of this derivative requires 5 x coordinates. One specific case is: Uxxxx(xg[i],yg[j]) = cxxxx[0]*U(xg[i-2],yg[j]) + cxxxx[1]*U(xg[i-1],yg[j]) + cxxxx[2]*U(xg[i ],yg[j]) + cxxxx[3]*U(xg[i+1],yg[j]) + cxxxx[4]*U(xg[i+2],yg[j]) The discrete derivative coefficients cxxxx[k] are computed based on xg[i-2]..xg[i+2] by a standard algorithm, implemented in many languages as nuderiv. The various U may be bound, known value, or free, unknown value, numeric values or variables. The variables will be determined by solving a system of linear simultaneous equations. Because a minimum of five vertices is required for fourth order discrete derivatives, a 5 by 5 array of vertices is the minimum that can be used for discretization of this example. Note that any of the 25 vertices may be bound or free. The array subscripts are chosen to start a zero. xg[i+0] xg[i+1] xg[i+2] xg[i+3] xg]i+4] +---------+---------+---------+---------+----------+ yg[j+0] | U(i,j) | | | | | +---------+---------+---------+---------+----------+ yg[j+1] | | | | | | +---------+---------+---------+---------+----------+ yg[j+2] | | | | | | +---------+---------+---------+---------+----------+ yg[j+3] | | | | | | +---------+---------+---------+---------+----------+ yg[j+4] | | | | |U(i+4,j+4)| +---------+---------+---------+---------+----------+ U(i,j) representing U(xg[i],yg[j]) Note that there are 25 possible locations for U(i,j) in the array. The offset on xg and yg change correspondingly: xg[i-2] xg[i-1] xg[i ] xg[i+1] xg]i+2] +---------+---------+---------+---------+----------+ yg[j-3] | | | | | | +---------+---------+---------+---------+----------+ yg[j-2] | | | | | | +---------+---------+---------+---------+----------+ yg[j-1] | | | | | | +---------+---------+---------+---------+----------+ yg[j ] | | | U(i,j) | | | +---------+---------+---------+---------+----------+ yg[j+1] | | | | | | +---------+---------+---------+---------+----------+ U(i,j) representing U(xg[i],yg[j]) The linear system of simultaneous equations that must be constructed is denoted A U = Y. A choice is possible, U may contain all vertices nx times ny values, or U may contain only the free vertices, DOF. The indexing for the latter choice is more complex, yet fewer equations need to be solved. There seems to be some, not large, difference in computation time and accuracy between the two choices. The A matrix will be nx*ny by nx*ny. The Y vector is the right hand side of the equations, f(xg[i],yg[j]) For each (i,j) vertex we will build an equation. The subscript i*ny+j is the equation number and row index. The column index will be computed from the subscripts of U. The full discrete equation for (i,j) in the (0,0) position of the 5 by 5 matrix, almost never written, is shown below. There are, of course, 25 such equations. Computer programs that implement discretization use variables and loops rather than implement the 25 equations. The coefficient arrays cxxxx, cyyyy, cxx, cyy and coefficient matrix cxxyy are computed based on xg[i]..xg[i+4] and yg[j]..yg[j+4]. The A[row][col] is summed with: row i*ny+j column (note pattern) cxxxx[0]* U(i ,j ) + (i )*ny+j summed with cxxxx[0] cxxxx[1]* U(i+1,j ) + (i+1)*ny+j summed with cxxxx[1] cxxxx[2]* U(i+2,j ) + cxxxx[3]* U(i+3,j ) + cxxxx[4]* U(i+4,j ) + 2*cxxyy[0][0]*U(i ,j ) + 2*cxxyy[0][1]*U(i+1,j ) + 2*cxxyy[0][2]*U(i+2,j ) + 2*cxxyy[0][3]*U(i+3,j ) + 2*cxxyy[0][4]*U(i+4,j ) + 2*cxxyy[1][0]*U(i ,j+1) + 2*cxxyy[1][1]*U(i+1,j+1) + (i+1)*ny+j+1 summed with 2*cxxyy[1][1] 2*cxxyy[1][2]*U(i+2,j+1) + (i+2)*ny+j+1 summed with 2*cxxyy[1][2] 2*cxxyy[1][3]*U(i+3,j+1) + 2*cxxyy[1][4]*U(i+4,j+1) + 2*cxxyy[2][0]*U(i ,j+2) + 2*cxxyy[2][1]*U(i+1,j+2) + 2*cxxyy[2][2]*U(i+2,j+2) + 2*cxxyy[2][3]*U(i+3,j+2) + 2*cxxyy[2][4]*U(i+4,j+2) + 2*cxxyy[3][0]*U(i ,j+3) + 2*cxxyy[3][1]*U(i+1,j+3) + 2*cxxyy[3][2]*U(i+2,j+3) + 2*cxxyy[3][3]*U(i+3,j+3) + 2*cxxyy[3][4]*U(i+4,j+3) + 2*cxxyy[4][0]*U(i ,j+4) + 2*cxxyy[4][1]*U(i+1,j+4) + 2*cxxyy[4][2]*U(i+2,j+4) + 2*cxxyy[4][3]*U(i+3,j+4) + 2*cxxyy[4][4]*U(i+4,j+4) + (i+4)*ny+j+4 cyyyy[0]* U(i ,j ) + (i )*ny+j cyyyy[1]* U(i ,j+1) + (i )*ny+j+1 cyyyy[2]* U(i ,j+2) + cyyyy[3]* U(i ,j+3) + cyyyy[4]* U(i ,j+4) + cxx [0]* U(i ,j ) + (i )*ny+j cxx [1]* U(i+1,j ) + (i+1)*ny+j cxx [2]* U(i+2,j ) + cxx [3]* U(i+3,j ) + cxx [4]* U(i+4,j ) + cyy [0]* U(i ,j ) + (i )*ny+j cyy [1]* U(i ,j+1) + (i )*ny+j+1 cyy [2]* U(i ,j+2) + cyy [3]* U(i ,j+3) + cyy [4]* U(i ,j+4) + 2 * U(i ,j ) (i )*ny+j summed with 2 This is repeated for all rows of the A matrix. Y[i*ny+j] is f(xg[i],yg[j]) The solution of the simultaneous equations computes U[i*ny+j] For the second choice, having equations only for free vertices, if the U(ik,jk) is a boundary, then the value of U(ik,jk) is known. Thus, subtract the coefficient times U(ik,jk) from Y, the right hand side. No action is taken for the A matrix for this case. Programming is tedious and requires attention to detail. There is no room for error. Any single subscript that is programmed incorrectly will cause incorrect results.