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);