From: "William R. Frensley" Subject: Re: Nine-point finite difference Laplace for irregular grid? Date: Wed, 15 Mar 2000 14:49:16 -0600 Newsgroups: sci.math.num-analysis Summary: [missing] moucha@helios.physics.utoronto.ca wrote: > > Hi, > I'm having trouble deriving a nine-point finite difference > representation for Laplace equation on an irregularly spaced > rectangular grid. > > I can derive the five-point for these types of grids, and a nine-point > for a regular square grid, using Taylor expansion. > Can someone point me in the right direction? Think of the problem like this: Your nine points are exact values for your function, and you want to fit a polynomial function to them. You can get a nine-term polynomial by taking all products of x^0, x^1, x^2 with y^0, y^1, y^2. When you have solved for the coefficients of this polynomial (probably a job for a symbolic algebra program), apply the continuum Laplace operator to the the polynomial and evaluate the result at the center point. This will give you the nine-point formula you need. A general hint: EVERY finite-difference formulation is based upon some assumptions about how the field varies between the mesh points. The thrust of the usual treatments of numerical methods is to try to ignore such details. Knowing exactly what the assumed form is will always pay off in the long run. - Bill Frensley ============================================================================== From: "James Van Buskirk" Subject: Re: Fixed nodes in Gauss quadratures Date: Sat, 7 Oct 2000 19:34:00 -0600 Newsgroups: sci.math.num-analysis Summary: [missing] Christian Proulx wrote in message <39DDFF20.6C3A0668@gmc.ulaval.ca>... >A classical Gauss quadrature is defined over ]-1,1[. The Gauss-Lobatto >quadrature is defined over [-1,1]. Is there a way to obtained >quadratures which would include to fixed absissas say -0.8 and 0.8? >I used the Fortran program DGQRUL, but I can only fix points outside the >[-1,1] interval. Why is that? If you think about what you are asking for, you want to be divide a polynomial of order at most 2*n+1 by P_n(x), which is the polynomial that has zeros at the n unknown points of the resulting integration formula: P_n(z) = product(/(z-x(i),i=1,n)/) times (x-a)*(x-b), where a and b are the two fixed points to get y_{2*n+1}(x) = Q_{n-1}(x)*P_n(x)*(x-a)*(x-b)+R_{n+1}(x) where Q_{n-1}(x) is the quotient of degree at most n-1 and R_{n+1}(x) is the remainder of degree at most n+1. For this to work, P_n(x) must be orthogonal to all polynomials of degree at most n-1 over the interval [-1,1] with weight function w(x) = (x-a)*(x-b) and have all of its roots simple and lying in the interval [-1,1]. This is not in general possible if either a or b lies in the interval (-1,1) because the weight function then no longer has the same sign throughout the interval. If you choose the weight function (x-a)**2 and use the value of y(a) and y'(a) in your integration formula then you can always get a valid Gaussian integration formula. Also for some values of a, b, and n you may still get a P_n(x) that has only simple roots lying in the interval [-1,1] and is orthogonal to all polynomials of degree at most n-1 over the interval [-1,1] with weight function w(x) = (x-a)*(x-b), even though a and b lie in (-1,1). Try a = -1+epsilon_1; b = 1-epsilon_2, where epsilon_1 and epsilon_2 are small and n not very large, either.