MATLAB三次样条插值之三弯矩法

ＭATLＡＢ三次样条插值之三弯矩法

首先说这个程序并不完善，为了实现通用（１，２，…,ｎ）格式解题，以及为调用追赶法程序，没有针对节点数在三个以下的情况进行分类讨论。希望能有朋友给出更好的方法。

首先，通过函数saｎwanｊ得到方程的系数矩阵，即追赶法方程的四个向量参数,接下来调

用追赶法(在inteｒｓaｎwj函数中），得到三次样条分段函数系数因子,然后进行多项式合并

得到分段函数的解析式，程序最后部分通过判断输入值的区间自动选择对应的分段函数并计算

改点的值。附:追赶法程序cｈase

％%%％％%%%%%%%％%

fｕncｔiｏｎ [nｅwv,w,newu，neｗd]=sanwj（ｘ,y,x０,y0,y1ａ,y1b）?％三弯矩样

条插值?％将插值点分两次输入，x0y０单独输入?% 边值条件a的二阶导数 y1a 和ｂ

的二阶导数ｙ1b，这里建议将y１ａ和y1b换成y2a和ｙ2b,以便于和三转角代码相区别

?n=length(x);m=leｎgth(y);

ｉf m~=n?errｏｒ(＇x or y 输入有误,再来'）;

end?v=ｏneｓ(ｎ-１，1);u=ｏnes(n－１,1）;d=zeros（n－１,１）;?w=２*ｏ

nes(n+１)；?h0=x（１）-x0;?h=zerｏs(n-１,1);

for k＝1:n-１?ｈ(k)＝x（k+1)－x（k);?enｄ

ｖ(1)＝h0/(ｈ0+h(1）);

u(1)＝1－v（1）;

ｄ（１)＝6*((y（2）-y(1）)/h（1)-(y(1)-y０)/ｈ0）/(h0+ｈ(１）);?%

ｆoｒ k＝2:ｎ-1?ｖ(ｋ）＝h(ｋ-１)/(ｈ(k-１)+h(ｋ)）;?u(k)=1－v(k);?d（k)＝

6*((y(ｋ+１)-ｙ（ｋ)）/h（k）-(ｙ(k）-y(ｋ-1）)/ｈ(k－1))／(h（k-1)+h(k));

ｅnd

newv=[v;1］;?newu=［1;u]；

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

d(n)=6＊（y1ｂ-(ｙ(n）－ｙ(n-1))／ｈ(n-1))／h(ｎ-1）;

newd=［d0;d];

%%%%%%％％%％％％

ｆunctｉon iｎｔersａnｗj（x，y,ｘ0,y0,y1a，y1b)

%三弯矩样条插值?%第一部分?n=leｎgtｈ(ｘ）;ｍ＝length(y);

iｆ m~＝n?eｒror（'ｘｏｒy 输入有误，再来');

enｄ?%重新定义h?h=zｅｒｏs(n,1）;

h(1)=x（1）-ｘ０;

fｏr k=2:n

h（k)=ｘ(ｋ)－x(k-1)；?end

％sptep1调用三弯矩函数?[a,b,c,d]=sanwj(x，ｙ,ｘ0,ｙ0,y１ａ,y1b);

% 三对角方程?M＝cｈaｓe(a,b,ｃ,ｄ);

% 求插值函数?fｐrintf（＇三次样条(三弯矩)插值的函数表达式\ｎ')；

sｙｍｓ X ；?fpｒintf('Ｓ0--1:\n');

S（１)＝cｏlｌeｃｔ（（(１/６）*M(2)*（X－x0)．^3－(1/6)*M(1)*(X－ｘ（1))．^3+(y(1)-(M(2)*h(１）.^2）/6)*(X－x０)－(y０－(Ｍ(1)*h（1)．^2）/6）＊(Ｘ-x(1））)/ｈ（1));?foｒ k=2：n

ｆpｒinｔｆ（＇Ｓ%d-－%d:\n',k-1,k)；

Ｓ(ｋ)=ｃｏllect((（1／６)*M(ｋ+1)*(X-x(k－1)).＾３-(1／6)*M(ｋ)*(X-x(ｋ)).^３＋(y(k）-（M(k+1)*h(k).^2)/6）*(X－ｘ（ｋ-1)）-(y(k-１)-(M(k)＊ｈ(k)．^2)/６)*(X-ｘ（k)))/h(ｋ)）;

ｅnd?S=S．';?ｄisp（S);?fprintf('以上为样条函数(三弯矩)解析式，显示为手写如下：＼n＇）；?preｔtｙ(S);

%第二部分

%是否继续运行程序?mｙlｏop=iｎput('继续运行程序输入“1”，否则输入“0”\n＇);?ｉｆmyloｏp

ｗhile ｍylｏop

xi＝input('输入需要计算的点的值，并按回车键\ｎ＇)；

if ｘi＞ｘ0｜xi

else?fprintｆ（'输入数值不在插值范围内,请重新输入\n');?ｘi=inpｕｔ('输入需要计算的点的值,并按回车键……＼ｎ');

eｎd

％确定输入的数值应该使用哪个解析式

newx=[x0;ｘ］;

[r,suoy］＝min(ａbs(newｘ－xi)）；?fpｒintｆ（'输入点的插值函数值为:\n\n＇)；fpr ｉntf(＇＼t＇);

if ｘi＜=nｅwx(sｕｏy）?f=ｓubs（S(suoy-１),X,xi）;?else

f=sｕｂs(S(suoy）,X,ｘｉ)；

eｎd

ｄisp(f);?mｙloop＝inｐut('继续计算输入“1”，终止计算输入“０”\n')；

end

else

return；

ｅnｄ

%%%%%%％%%%%%

functｉon[x]=chasｅ(a,ｂ,c,d)

%追赶法解性方程组 a是下三角b是对角线c是上三角 d是常数项

%输入的a ｂｃ d 均为列向量?n=lengｔh(b)；?u=ｚｅros(n,1）;?v=zeroｓ(n，1);?x=zerｏs（n,1);

%追

v(１）=c(1）/b(1);u(１)=d(1)/b(１）；

for i＝2：n-１

v（i)=c（ｉ)／（b(i)-v(i－1)*a(ｉ－1)）；

u(i)=(d(i)－u（i-１)*a（i-1))/(b(i)-v(i-１)＊a(i-1）)；

enｄ?u(n）＝(ｄ(n）－ｕ(n-1）*a(n－1))/（b（n)-v(n-１)*a(n-1）)；?%赶?x （n)=u(n）;?foｒｉ=n-1:-１:1

x(i)＝u(i)－v(ｉ)*ｘ(ｉ+1);?eｎd