From: "William R. Frensley" <frensley@utdallas.edu>
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" <torsop@ix.netcom.com>
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.