diff options
author | jthorn <jthorn@df1f8a13-aa1d-4dd4-9681-27ded5b42416> | 2002-08-28 11:05:53 +0000 |
---|---|---|
committer | jthorn <jthorn@df1f8a13-aa1d-4dd4-9681-27ded5b42416> | 2002-08-28 11:05:53 +0000 |
commit | afcc4e7c56586ac436b81763ab145d1ef3feea65 (patch) | |
tree | f61ef4ba6ebc77e4bd76a95cbf5efe794b4d98e5 /src/GeneralizedPolynomial-Uniform/interpolate.maple | |
parent | 354d16ee3d9033f468474fcd4b71cffeccf947ae (diff) |
redo Hermite interpolant computation to hopefully be more efficient
(in practice it doesn't seem to make much difference, though :( :(),
and to be a bit easier to read/understand/modify the code
git-svn-id: http://svn.cactuscode.org/arrangements/CactusBase/LocalInterp/trunk@99 df1f8a13-aa1d-4dd4-9681-27ded5b42416
Diffstat (limited to 'src/GeneralizedPolynomial-Uniform/interpolate.maple')
-rw-r--r-- | src/GeneralizedPolynomial-Uniform/interpolate.maple | 198 |
1 files changed, 105 insertions, 93 deletions
diff --git a/src/GeneralizedPolynomial-Uniform/interpolate.maple b/src/GeneralizedPolynomial-Uniform/interpolate.maple index 3f627db..58a2c97 100644 --- a/src/GeneralizedPolynomial-Uniform/interpolate.maple +++ b/src/GeneralizedPolynomial-Uniform/interpolate.maple @@ -86,7 +86,7 @@ end proc; # # This function computes a Hermite polynomial interpolant in any # number of dimensions. This is a polynomial which -# has values which match the given data DATA() at a specified set of +# * 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 @@ -105,50 +105,45 @@ end proc; # + 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 +# coeff_set = A set of the interpolation coefficients (coefficients in # 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]. -# 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 +# deriv_set = A set of equations of the form +# {coords} = proc +# giving the derivatives which are to be matched, and the +# procedures to compute their finite-difference approximations. +# 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 +# { +# {x} = 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 +# {y} = 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]] +# {x,y} = 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_set = A set 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_set = A list of positions (each a list of numeric values) +# where the interpolant is to match the derivatives +# specified by deriv_set , for example +# {[0,0], [0,1], [1,0], [1,1]} # # Results: # This function returns the interpolating polynomial, in the form of @@ -156,71 +151,88 @@ end proc; # Hermite_polynomial_interpolant := proc( - 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)) + fn::procedure, + coeff_set::set(name), + coord_list::list(name), + deriv_set::set(set(name) = procedure), + fn_posn_set::set(list(numeric)), + deriv_posn_set::set(list(numeric)) ) -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)); +local fn_eqnset, deriv_eqnset, coeff_eqns, subs_eqnset; + + +# +# compute a set of equations +# {fn(posn) = DATA(posn)} +# giving the function values to be matched +# +fn_eqnset := map( + # return equation that fn(posn) = DATA(posn) + proc(posn::list(integer)) + fn(op(posn)) = 'DATA'(op(posn)); + end proc + , + fn_posn_set + ); + + +# +# compute a set of equations +# { diff(fn,coords)(posn) = DERIV(coords)(posn) } +# giving the derivative values to be matched, where DERIV(coords) +# is a placeholder for the appropriate derivative +# 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)] - ) - ) - } + # return set of equations for this particular derivative + proc(deriv_coords::set(name)) + local deriv_fn; + fn(op(coord_list)); + diff(%, op(deriv_coords)); + deriv_fn := unapply(%, op(coord_list)); + map( + proc(posn::list(integer)) + deriv_fn(op(posn)) = 'DERIV'(op(deriv_coords))(op(posn)); + end proc + , + deriv_posn_set + ); end proc , - {op(deriv_posn_list)} + map(lhs, deriv_set) ); +deriv_eqnset := `union`(op(%)); + + +# +# solve overall set of equations for coefficients +# in terms of DATA() and DERIV() values +# +coeff_eqns := solve[linear](fn_eqnset union deriv_eqnset, coeff_set); + + +# +# compute a set of substitution equations +# {'DERIV'(coords) = procedure} +# +subs_eqnset := map( + proc(eqn::set(name) = procedure) + 'DERIV'(op(lhs(eqn))) = rhs(eqn); + end proc + , + deriv_set + ); -# 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)}); +# +# compute the coefficients in terms of the DATA() values +# +subs(subs_eqnset, coeff_eqns); +eval(%); -# interpolant as a polynomial in the coordinates -return subs(%, eval(fn))(op(coord_list)); +# +# compute the interpolant as a polynomial in the coordinates +# +subs(%, fn(op(coord_list))); end proc; ################################################################################ |