我正在 matlab 中编写一个算法来对表面进行双三次插值Psi(x,y)。我在代码中有一个错误，似乎无法追踪它。我正在尝试一个测试用例，Psi=X^2-0.25以便更容易追踪错误。好像我的插值有一个偏移量。我的评论包含在我的代码中。任何帮助，将不胜感激。
蓝色的绘图Psi=X^2和红色的插值
绘制了Countour 线，Psi红点是我正在计算插值的点。粗红线是从红点偏移相当多的插值。
function main()
epsilon=0.000001;
xMin=-1+epsilon;
xMax= 1+epsilon;
yMin=-1+epsilon;
yMax= 1+epsilon;
dx=0.1; Nx=ceil((xMax-xMin)/dx)+1;
dy=0.1; Ny=ceil((yMax-yMin)/dy)+1;
x=xMin:dx:xMax; x=x(1:Nx);
y=yMin:dy:yMax; y=y(1:Ny);
[XPolInX,XPolInY]=GetGhostMatricies(Nx,Ny); %Linear extrapolation matrix
[D0x,D0y]=GetDiffMatricies(Nx,Ny,dx,dy); %derivative matricies: D0x is central differencing in x
[X,Y]=meshgrid(x,y);
Psi=X.^2-0.25; %Note that my algorithm is being written for a Psi that may not have a analytic representation. This Psi is only a test case.
psi=zeros(Nx+2,Ny+2); %linearly extrapolate psi (for solving differential equation not shown here)
psi(2:(Nx+1),2:(Ny+1))=Psi';
psi=(XPolInY*(XPolInX*psi)')';
%compute derivatives of psi
psi_x =D0x*psi; psi_x =(XPolInY*(XPolInX*psi_x)')';
psi_y =(D0y*psi')'; psi_y =(XPolInY*(XPolInX*psi_y)')';
psi_xy=D0x*psi_y; psi_xy=(XPolInY*(XPolInX*psi_xy)')';
% i have verified that my derivatives are computed correctly
biCubInv=GetBiCubicInverse(dx,dy);
i=5; %lets compute the bicubic interpolation at this x(i), y(j)
j=1;
psiVoxel=[psi( i,j),psi( i+1,j),psi( i,j+1),psi( i+1,j+1),...
psi_x( i,j),psi_x( i+1,j),psi_x( i,j+1),psi_x( i+1,j+1),...
psi_y( i,j),psi_y( i+1,j),psi_y( i,j+1),psi_y( i+1,j+1),...
psi_xy(i,j),psi_xy(i+1,j),psi_xy(i,j+1),psi_xy(i+1,j+1)]';
a=biCubInv*psiVoxel; %a=[a00 a01 ... a33]; polynomial coefficients; 1st index is power of (x-xi), 2nd index is power of (y-yj)
xi=x(5); yj=y(1);
clear x y
x=(xi-.2):.01:(xi+.2); %this is a local region about the point we are interpolating
y=(yj-.2):.01:(yj+.2);
[dX,dY]=meshgrid(x,y);
Psi=dX.^2-0.25;
figure(2) %just plotting the 0 level contour of Psi here
plot(xi,yj,'.r','MarkerSize',20)
hold on
contour(x,y,Psi,[0 0],'r','LineWidth',2)
set(gca,'FontSize',14)
axis([x(1) x(end) y(1) y(end)])
grid on
set(gca,'xtick',(xi-.2):.1:(xi+.2));
set(gca,'ytick',(yj-.2):.1:(yj+.2));
xlabel('x')
ylabel('y')
[dX dY]=meshgrid(x-xi,y-yj);
%P is my interpolating polynomial
P = a(1) + a(5) *dY + a(9) *dY.^2 + a(13) *dY.^3 ...
+ a(2)*dX + a(6)*dX .*dY + a(10)*dX .*dY.^2 + a(14)*dX .*dY.^3 ...
+ a(3)*dX.^2 + a(7)*dX.^2.*dY + a(11)*dX.^2.*dY.^2 + a(15)*dX.^2.*dY.^3 ...
+ a(4)*dX.^3 + a(8)*dX.^3.*dY + a(12)*dX.^3.*dY.^2 + a(16)*dX.^3.*dY.^3 ;
[c h]=contour(x,y,P)
clabel(c,h)
figure(3)
plot(x,x.^2-.25) %this is the exact function
hold on
plot(x,P(1,:),'-r*')
%See there is some offset here
end
%-------------------------------------------------------------------------
function [XPolInX,XPolInY]=GetGhostMatricies(Nx,Ny)
XPolInX=diag(ones(1,Nx+2),0);
XPolInY=diag(ones(1,Ny+2),0);
XPolInX(1,1) =0; XPolInX(1,2) =2; XPolInX(1,3) =-1;
XPolInY(1,1) =0; XPolInY(1,2) =2; XPolInY(1,3) =-1;
XPolInX(Nx+2,Nx+2)=0; XPolInX(Nx+2,Nx+1)=2; XPolInX(Nx+2,Nx)=-1;
XPolInY(Ny+2,Ny+2)=0; XPolInY(Ny+2,Ny+1)=2; XPolInY(Ny+2,Ny)=-1;
fprintf('Done GetGhostMatricies\n')
end
%-------------------------------------------------------------------------
function [D0x,D0y]=GetDiffMatricies(Nx,Ny,dx,dy)
D0x=diag(ones(1,Nx-1),1)-diag(ones(1,Nx-1),-1);
D0y=diag(ones(1,Ny-1),1)-diag(ones(1,Ny-1),-1);
D0x(1,1)=-3; D0x(1,2)=4; D0x(1,3)=-1;
D0y(1,1)=-3; D0y(1,2)=4; D0y(1,3)=-1;
D0x(Nx,Nx)=3; D0x(Nx,Nx-1)=-4; D0x(Nx,Nx-2)=1;
D0y(Ny,Ny)=3; D0y(Ny,Ny-1)=-4; D0y(Ny,Ny-2)=1;
%pad with ghost cells which are simply zeros
tmp=D0x; D0x=zeros(Nx+2,Nx+2); D0x(2:(Nx+1),2:(Nx+1))=tmp; tmp=0;
tmp=D0y; D0y=zeros(Ny+2,Ny+2); D0y(2:(Ny+1),2:(Ny+1))=tmp; tmp=0;
%scale appropriatley by dx & dy
D0x=D0x/(2*dx);
D0y=D0y/(2*dy);
end
%-------------------------------------------------------------------------
function biCubInv=GetBiCubicInverse(dx,dy)
%p(x,y)=a00+a01(x-xi)+a02(x-xi)^2+...+a33(x-xi)^3(y-yj)^3
%biCubic*a=[psi(i,j) psi(i+1,j) psi(i,j+1) psi(i+1,j+1) psi_x(i,j) ... psi_y(i,j) ... psi_xy(i,j) ... psi_xy(i+1,j+1)]
%here, psi_x is the x derivative of psi
%I verified that this matrix is correct by setting dx=dy=1 and comparing to the inverse here http://en.wikipedia.org/wiki/Bicubic_interpolation
biCubic=[
%00 10 20 30 01 11 21 31 02 12 22 32 03 13 23 33
1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
1 dx dx^2 dx^3 0 0 0 0 0 0 0 0 0 0 0 0;
1 0 0 0 dy 0 0 0 dy^2 0 0 0 dy^3 0 0 0;
1 dx dx^2 dx^3 dy dx*dy dx^2*dy dx^3*dy dy^2 dx*dy^2 dx^2*dy^2 dx^3*dy^2 dy^3 dx*dy^3 dx^2*dy^3 dx^3*dy^3;
0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 2*dx 3*dx^2 0 0 0 0 0 0 0 0 0 0 0 0;
0 1 0 0 0 dy 0 0 0 dy^2 0 0 0 dy^3 0 0;
0 1 2*dx 3*dx^2 0 dy 2*dx*dy 3*dx^2*dy 0 dy^2 2*dx*dy^2 3*dx^2*dy^2 0 dy^3 2*dx*dy^3 3*dx^2*dy^3;
0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 1 dx dx^2 dx^3 0 0 0 0 0 0 0 0;
0 0 0 0 1 0 0 0 2*dy 0 0 0 3*dy^2 0 0 0;
0 0 0 0 1 dx dx^2 dx^3 2*dy 2*dx*dy 2*dx^2*dy 2*dx^3*dy 3*dy^2 3*dx*dy^2 3*dx^2*dy^2 3*dx^3*dy^2;
0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 2*dx 3*dx^2 0 0 0 0 0 0 0 0;
0 0 0 0 0 1 0 0 0 2*dy 0 0 0 3*dy^2 0 0;
0 0 0 0 0 1 2*dx 3*dx^2 0 2*dy 4*dx*dy 6*dx^2*dy 0 3*dy^2 6*dx*dy^2 9*dx^2*dy^2];
biCubInv=inv(biCubic);
end
%-------------------------------------------------------------------------