switch from $Id:$ to $Header:$

git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@93 df1f8a13-aa1d-4dd4-9681-27ded5b42416

1 files changed, 119 insertions, 20 deletions

diff --git a/src/GeneralizedPolynomial-Uniform/interpolate.maple b/src/GeneralizedPolynomial-Uniform/interpolate.maple
index 23ca9e0..3f627db 100644
--- a/src/GeneralizedPolynomial-Uniform/interpolate.maple
+++ b/src/GeneralizedPolynomial-Uniform/interpolate.maple
@@ -1,5 +1,5 @@
 # interpolate.maple -- compute interpolation formulas/coefficients
-# $Id$
+# $Header$
 
 #
 # <<<representation of numbers, data values, etc>>>
@@ -85,21 +85,70 @@ end proc;
 
 #
 # This function computes a Hermite polynomial interpolant in any
-# number of dimensions.
+# number of dimensions.  This is a polynomial which
+#  has values which match the given data DATA() at a specified set of
+#   points, and
+# * has derivatives which match the specified finite-difference derivatives
+#   of the given data DATA() at a specified set of points
+#
+# For the derivative matching, we actually match all possible products
+# of 1st derivatives, i.e. in 2-D we match dx, dy, and dxy, in 3-D we
+# match dx, dy, dz, dxy, dxz, dyz, and dxyz, etc etc.
 #
 # Arguments:
 # fn = The interpolation function.  This should be a procedure in the
 #      coordinates, having the coefficients as global variables.  For
 #      example,
-#	  proc(x,y) c00 + c10*x + c01*y end proc
+#		proc(x,y)
+#		+ c03*y^3 + c13*x*y^3 + c23*x^2*y^3 + c33*x^3*y^3
+#		+ c02*y^2 + c12*x*y^2 + c22*x^2*y^2 + c32*x^3*y^2
+#		+ c01*y   + c11*x*y   + c21*x^2*y   + c31*x^3*y
+#		+ c00     + c10*x     + c20*x^2     + c30*x^3
+#		end proc;
 # coeff_list = A set of the interpolation coefficients (coefficients in
-#	       the interpolation function), for example [c00, c10, c01].
+#	       the interpolation function), for example
+#			[
+#			c03, c13, c23, c33,
+#			c02, c12, c22, c32,
+#			c01, c11, c21, c31,
+#			c00, c10, c20, c30
+#			]
 # coord_list = A list of the coordinates (independent variables in the
 #	       interpolation function), for example [x,y].
-# posn_list = A list of positions (each a list of numeric values) where the
-#	      interpolant is to use data, for example  hypercube([0,0], [1,1]).
-#	      Any positions may be used; if they're redundant (as in the
-#	      example) the least-squares interpolant is computed.
+# deriv_coord_list = A list of lists of coordinates specifying which
+#		     derivatives are computed by the deriv_proc_list[]
+#		     procedures, for example
+#			[[x], [y], [x,y]]
+# deriv_proc_list = A list of procedures for computing finite-difference
+#		    derivatives.  Each procedure should take N_dims integer
+#		    arguments specifying an evaluation point, and return
+#		    a suitable linear combination of the DATA() for the
+#		    derivative at that point.  For example
+#		    example,
+#			[
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i-1,j) + 1/2*DATA(i+1,j)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/2*DATA(i,j-1) + 1/2*DATA(i,j+1)
+#			  end proc
+#			,
+#			  proc(i::integer, j::integer)
+#			  - 1/4*DATA(i-1,j+1) + 1/4*DATA(i+1,j+1)
+#			  + 1/4*DATA(i-1,j-1) - 1/4*DATA(i+1,j-1)
+#			  end proc
+#			]
+# fn_posn_list = A list of positions (each a list of numeric values)
+#		 where the interpolant is to match the given data DATA(),
+#		 for example
+#			[[0,0], [0,1], [1,0], [1,1]]
+# deriv_posn_list = A list of positions (each a list of numeric values)
+#		    where the interpolant is to match each possible product
+#		    of 1st derivatives"1st"
+#		    derivatives (actually all derivatives the given data DATA(),
+#		    for example
+#			[[0,0], [0,1], [1,0], [1,1]]
 #
 # Results:
 # This function returns the interpolating polynomial, in the form of
@@ -107,20 +156,71 @@ end proc;
 #
 Hermite_polynomial_interpolant :=
 proc(
-      fn::procedure, coeff_list::list(name),
-      coord_list::list(name), posn_list::list(list(numeric))
+      fn::procedure, coeff_list::list(name), coord_list::list(name),
+      deriv_coord_list::list(list(name)), fd_deriv_proc_list::list(procedure),
+      fn_posn_list::list(list(numeric)), deriv_posn_list::list(list(numeric))
     )
-local posn, data_eqns, coeff_eqns;
+local fn_eset,
+      fn_expr, deriv_eset;
+
+# set of equations {fn(posn) = DATA(posn)}
+fn_eset := map(
+		  # return equation that fn(this point) = DATA(this point)
+		  proc(posn::list(integer))
+		  fn(op(posn)) = 'DATA'(op(posn));
+		  end proc
+		,
+		  {op(fn_posn_list)}
+	      );
+
+# set of sets of equations
+#	{
+#	  { deriv1(posn1) = fd_deriv1(posn1),
+#	    deriv2(posn1) = fd_deriv2(posn1),
+#	    ... },
+#	  { deriv1(posn2) = fd_deriv1(posn2),
+#	    deriv2(posn2) = fd_deriv2(posn2),
+#	    ... },
+#	  ...
+#	}
+fn_expr := fn(op(coord_list));
+map(
+       # return set of equations
+       #	{deriv(posn) = fd_deriv(posn)}
+       # for this point
+       proc(posn::list(integer))
+       {
+	op(
+	 zip(
+		# return equation that
+		#	deriv(posn) = fd_deriv(posn)
+		# for this deriv and this point
+		proc(deriv_coords::list(name), fd_deriv_proc::procedure)
+		local fn_deriv_proc;
+		fn_deriv_proc := unapply( diff(fn_expr,deriv_coords),
+					  op(coord_list) );
+		fn_deriv_proc(op(posn)) = fd_deriv_proc(op(posn));
+		end proc
+	      ,
+		[op(deriv_coord_list)]
+	      ,
+		[op(fd_deriv_proc_list)]
+	    )
+	  )
+       }
+       end proc
+     ,
+       {op(deriv_posn_list)}
+   );
 
-# coefficients of interpolating polynomial
-data_eqns := {  seq( fn(op(posn))='DATA'(op(posn)) , posn=posn_list )  };
-coeff_eqns := linalg[leastsqrs](data_eqns, {op(coeff_list)});
-if (has(coeff_eqns, '_t'))
-   then error "interpolation coefficients aren't uniquely determined!";
-end if;
+# set of equations {deriv-i(posn-j) = fd_deriv-i(posn-j)}
+deriv_eset := `union`(op(%));
+
+# solve equations for coeffs
+solve(fn_eset union deriv_eset, {op(coeff_list)});
 
 # interpolant as a polynomial in the coordinates
-return subs(coeff_eqns, eval(fn))(op(coord_list));
+return subs(%, eval(fn))(op(coord_list));
 end proc;
 
 ################################################################################
@@ -133,8 +233,7 @@ end proc;
 # Arguments:
 # interpolant = The interpolating polynomial (an algebraic expression
 #		in the coordinates and the data values).
-# posn_list = The same list of positions as was used to compute the
-#	      interpolating polynomial.
+# posn_list = The same list of data positions used in the interpolant.
 #
 # Results:
 # This function returns the coefficients, as a list of equations of the