Mikhail Genkin
Mikhail Genkin

Reputation: 3460

Return the coefficients of bicubic interpolation

Suppose, I have a function z=f(x,y). I calculated it in 16 points, say, at the square [0:3,0:3] (x=[0,1,2,3] y=[0,1,2,3]). Now I want to do bicubic interpolation. But as a result, I need the coefficients of this interpolation.

interp2 lets me to calculate the result of interpolation in the predefined collocation points, but it doesn't return the coefficients itself.

spline works only for 1D functions.

Any idea how I can get the coefficients of bicubic interpolation? I need the coefficients, because after that I intend to use them for the numerical solver (and I need the analytical expression of my function for that).

Upvotes: 0

Views: 716

Answers (1)

Mikhail Genkin
Mikhail Genkin

Reputation: 3460

I just wrote down the function that returns coefficients. Assuming 16 values comes from (x=-dx,0,dx,2dx) and (y=-dx,0,dx,2dx). I also wrote function 'Deriv' that takes derivatives of my numerical function with respect to 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

Upvotes: 1

Related Questions