MATLAB三次樣條插值之三彎矩法1

第頁MATLAB三次樣條插值之三彎矩法首先說這個(gè)程序并不完善，為了實(shí)現(xiàn)通用（1，2，…，n）格式解題，以及為調(diào)用追趕法程序，沒有針對節(jié)點(diǎn)數(shù)在三個(gè)以下的情況進(jìn)行分類討論。希望能有朋友給出更好的方法。首先，通過函數(shù)sanwanj得到方程的系數(shù)矩陣，即追趕法方程的四個(gè)向量參數(shù)，接下來調(diào)用追趕法（在intersanwj函數(shù)中），得到三次樣條分段函數(shù)系數(shù)因子，然后進(jìn)行多項(xiàng)式合并得到分段函數(shù)的解析式，程序最后部分通過判斷輸入值的區(qū)間自動(dòng)選擇對應(yīng)的分段函數(shù)并計(jì)算改點(diǎn)的值。附：追趕法程序chasefunction[newv,w,newu,newd]=sanwj(x,y,x0,y0,y1a,y1b)

%三彎矩樣條插值

%將插值點(diǎn)分兩次輸入，x0y0單獨(dú)輸入

%邊值條件a的二階導(dǎo)數(shù)y1a和b的二階導(dǎo)數(shù)y1b，這里建議將y1a和y1b換成y2a和y2b，以便于和三轉(zhuǎn)角代碼相區(qū)別

n=length(x);m=length(y);

ifm~=n

error('xory輸入有誤，再來');

end

v=ones(n-1,1);u=ones(n-1,1);d=zeros(n-1,1);

w=2*ones(n+1);

h0=x(1)-x0;

h=zeros(n-1,1);

fork=1:n-1

h(k)=x(k+1)-x(k);

end

v(1)=h0/(h0+h(1));

u(1)=1-v(1);

d(1)=6*((y(2)-y(1))/h(1)-(y(1)-y0)/h0)/(h0+h(1));

%

fork=2:n-1

v(k)=h(k-1)/(h(k-1)+h(k));

u(k)=1-v(k);

d(k)=6*((y(k+1)-y(k))/h(k)-(y(k)-y(k-1))/h(k-1))/(h(k-1)+h(k));

end

newv=[v;1];

newu=[1;u];

d0=6*((y(1)-y0)/h0-y1a)/h0;

d(n)=6*(y1b-(y(n)-y(n-1))/h(n-1))/h(n-1);

newd=[d0;d];functionintersanwj(x,y,x0,y0,y1a,y1b)

%三彎矩樣條插值

%第一部分

n=length(x);m=length(y);

ifm~=n

error('xory輸入有誤，再來');

end

%重新定義h

h=zeros(n,1);

h(1)=x(1)-x0;

fork=2:n

h(k)=x(k)-x(k-1);

end

%sptep1調(diào)用三彎矩函數(shù)

[a,b,c,d]=sanwj(x,y,x0,y0,y1a,y1b);

%三對角方程

M=chase(a,b,c,d);

%求插值函數(shù)

fprintf('三次樣條(三彎矩)插值的函數(shù)表達(dá)式\n');

symsX;

fprintf('S0--1:\n');

S(1)=collect(((1/6)*M(2)*(X-x0).^3-(1/6)*M(1)*(X-x(1)).^3+(y(1)-(M(2)*h(1).^2)/6)*(X-x0)-(y0-(M(1)*h(1).^2)/6)*(X-x(1)))/h(1));

fork=2:n

fprintf('S%d--%d:\n',k-1,k);

S(k)=collect(((1/6)*M(k+1)*(X-x(k-1)).^3-(1/6)*M(k)*(X-x(k)).^3+(y(k)-(M(k+1)*h(k).^2)/6)*(X-x(k-1))-(y(k-1)-(M(k)*h(k).^2)/6)*(X-x(k)))/h(k));

end

S=S.';

disp(S);

fprintf('以上為樣條函數(shù)(三彎矩)解析式，顯示為手寫如下：\n');

pretty(S);

%第二部分

%是否繼續(xù)運(yùn)行程序

myloop=input('繼續(xù)運(yùn)行程序輸入“1”，否則輸入“0”\n');

ifmyloop

whilemyloop

xi=input('輸入需要計(jì)算的點(diǎn)的值，并按回車鍵\n');

ifxi>x0|xi<x(n)

fprintf('現(xiàn)在開始計(jì)算輸入點(diǎn)的插值函數(shù)值……\n');

else

fprintf('輸入數(shù)值不在插值范圍內(nèi)，請重新輸入\n');

xi=input('輸入需要計(jì)算的點(diǎn)的值，并按回車鍵……\n');

end

%確定輸入的數(shù)值應(yīng)該使用哪個(gè)解析式

newx=[x0;x];

[r,suoy]=min(abs(newx-xi));

fprintf('輸入點(diǎn)的插值函數(shù)值為:\n\n');fprintf('\t');

ifxi<=newx(suoy)

f=subs(S(suoy-1),X,xi);

else

f=subs(S(suoy),X,xi);

end

disp(f);

myloop=input('繼續(xù)計(jì)算輸入“1”，終止計(jì)算輸入“0”\n');

end

else

return;

endfunction[x]=chase(a,b,c,d)

%追趕法解性方程組a是下三角b是對角線c是上三角d是常數(shù)項(xiàng)

%輸入的abcd均為列向量

n=length(b);

u=zeros(n,1);

v=zeros(n,1);

x=zeros(n,1);

%追

v(1)=c(1)/b(1);u(1)=d(1)/b(1);

fori=2:n-1

v(i)=c(i)/(b(i)-v(i-1)*a(i-1));

u(i)=(d(i)-u(i-1)*a(i-1))/(b(i)-v