and the inverse is
[9]
X
−1
=
1
4
−
1
2
1
4
−
1
2
1 −
1
2
1
4
−
1
2
1
4
−
3
4
1 −
1
4
3
2
−2
1
2
−
3
4
1 −
1
4
−
3
4
3
2
−
3
4
1 −21−
1
4
1
2
−
1
4
1
2
00−10 0
1
2
00
1
2
−1
1
2
000000
9
4
−3
3
4
−34−1
3
4
−1
1
4
−
3
2
00200−
1
2
00
−
3
2
2 −
1
2
000000
100000000
and the coefficient vector solution a is found simply by
[10]
a = X
−1
$ z
Again X
-1
is a constant matrix, so the coefficients for the interpolating polynomial are found with a single
matrix multiply. Another advantage to this method is that all of the elements of the inverted matrix can
be represented exactly with 89/92+ BCD arithmetic, so they do not contribute to round-off error in the
final result
.Once we have the polynomial coefficients, it is straight-forward to interpolate for z with
[11]
z = a $ u
2
$ v
2
+ b $ u
2
$ v + c $ u $ v
2
+ d $ u
2
+ e $ v
2
+ f $ u $ v + g $ u + h $ v + i
This function, intrp9z(), implements these ideas.
intrpz9(xl,yl,zmat,x,y)
Func
©({xlist},{ylist},[zmatrix],x,y) 9-point z-interpolation
©Uses matrix math\im1a
©1apr01/dburkett@infinet.com
local u,v
when(x<xl[2],(x-xl[1])/(xl[2]-xl[1]),(x-xl[2])/(xl[3]-xl[2])+1)→u
when(y<yl[2],(y-yl[1])/(yl[2]-yl[1]),(y-yl[2])/(yl[3]-yl[2])+1)→v
sum(mat▶list(math\im1a*(augment(augment(zmat[1],zmat[2]),zmat[3])))*{u^2*v^2,u^
2*v,u*v^2,u^2,v^2,u*v,u,v,1})
EndFunc
Note that the matrix im1a (from equation [9]) must be present and stored in the \math folder.
The input arguments are
6 - 50