2017-02-25 77 views
0

假設我有一個函數z=f(x,y)。我在16點計算它,例如在方塊[0:3,0:3]x=[0,1,2,3]y=[0,1,2,3])。現在我想做雙三次插值。但結果是,我需要這種插值的係數。返回雙三次插值係數

interp2讓我來計算預定義搭配點的插值結果,但它不返回係數本身。

spline僅適用於1D功能。

任何想法如何我可以得到雙三次插值的係數? 我需要係數,因爲之後我打算將它們用於數值解算器(並且我需要我的函數的解析表達式)。

回答

0

我剛寫下了返回係數的函數。假設16個值來自(x = -dx,0,dx,2dx)和(y = -dx,0,dx,2dx)。我還寫了函數'Deriv',它涉及我的數值函數關於x,y,xy的導數:

function [out]=Deriv(in,str,dx) 
    switch str 
    case 'x' 
     out=(circshift(in,[0 -1])-circshift(in,[0 1]))./(2*dx); 
    case 'y'  
     out=(circshift(in,[-1 0])-circshift(in,[1 0]))./(2*dx); 
    case 'xy' 
     out=(circshift(in,[-1 -1])+circshift(in,[1 1])-circshift(in,[-1 1])-circshift(in,[1 -1]))./(4*dx*dx); 
    end 

function out=bl(in1,in2) 
    out=[in1^3, in1^2, in1, 1,in2*[in1^3, in1^2, in1, 1], in2^2*[in1^3, in1^2, in1, 1]... 
    in2^3*[in1^3, in1^2, in1, 1]]; 
end 

function out=blx(in1,in2) 
    out=[3*in1^2, 2*in1, 1, 0,in2*[3*in1^2, 2*in1, 1, 0], in2^2*[3*in1^2, 2*in1, 1, 0]... 
    in2^3*[3*in1^2, 2*in1, 1, 0]]; 
end 

function out=bly(in1,in2) 
    out=[0,0,0,0,[in1^3, in1^2, in1, 1], 2*in2*[in1^3, in1^2, in1, 1]... 
    3*in2^2*[in1^3, in1^2, in1, 1]]; 
end 


function out=MyBicubic(in,dx) 
A=zeros(16,16); C=zeros(16,1); 
    A(1,:)=bl(0,0); C(1)=in(2,2); 
    A(2,:)=bl(dx,0); C(2)=in(2,3); 
    A(3,:)=bl(0,dx); C(3)=in(3,2); 
    A(4,:)=bl(dx,dx); C(4)=in(3,3); 
    A(5,:)=blx(0,0); temp=Deriv(in,'x',dx); C(5)=temp(2,2); 
    A(6,:)=blx(dx,0); C(6)=temp(2,3); 
    A(7,:)=blx(0,dx); C(7)=temp(3,2); 
    A(8,:)=blx(dx,dx); C(8)=temp(3,3); 
    A(9,:)=bly(0,0); temp=Deriv(in,'y',dx); C(9)=temp(2,2); 
    A(10,:)=bly(dx,0); C(10)=temp(2,3); 
    A(11,:)=bly(0,dx); C(11)=temp(3,2); 
    A(12,:)=bly(dx,dx); C(12)=temp(3,3); 
    A(13,:)=blxy(0,0); temp=Deriv(in,'xy',dx); C(13)=temp(2,2); 
    A(14,:)=blxy(dx,0); C(14)=temp(2,3); 
    A(15,:)=blxy(0,dx); C(15)=temp(3,2); 
    A(16,:)=blxy(dx,dx); C(16)=temp(3,3); 
    out=A\C; 
end