三次样条插值

三次样条插值1.算法原理由于在许多实际问题中，要求函数的二阶导数连续，人们便提出了三次样条插值函数，三次样条插值函数是由分段三次函数拼接而成的，在连接点处二阶导数连续。设S(x)在节点ix处的二阶导数),1,0()(''niMxSii，，其中iM为待定参数。由S（x）是分段三次多项式可知，)(''xS是分段线性函数，)(''xS在子区间iixx,1上可以表示为iiiiiiiiiiiiiiiixxxMhxxMhxxMxxxxMxxxxxS1111111,)(''其中1iiixxh，对S（x）两端积分两次得iiiiiiiiiiiixxxxxcxxbMhxxMhxxxS111113),(6)(6)()(其中ib和ic为积分常数。由插值条件iiiiyxSyxS)(,11得iiiiiiiiiiyhcMhyhbMh6,62112由此解得iiiiiiiiiihMhychMhyb/6,/62121代入得：iiiiiiiiiiiiiiiiiixxxhxxMhyhxxMhyMhxxMhxxxS1121213113,6666求导得：iiiiiiiiiiiiiixxxhMMhyyMhxxMhxxxS1112112,622'令ixx得xS在ix处的左导数iiiiiiiiiiiiiiihyyMhMhhMMhyyMhxS11113662'又令1ixx得xS在1ix处右导数iiiiiiiiiiiiiiiihyyMhMhhMMhyyMhxS1111116362'，从而有1111163)('iiiiiiihyyMhMhxS，由xS在节点ix处一阶导数的连续性知1,2,1),(')('nixSxSii，，即1,2,16361111111nihyyhyyMhMhhMhiiiiiiiiiiiii，，两端同乘16iihh得)(62111111111iiiiiiiiiiiiiiiiihyyhyyhhMhhhMMhhh，记iiiiiiiiihhhhhh1,111，1,,2,1],,,[66111111nixxxfhyyhyyhhdiiiiiiiiiiii，则关于iM的方程组写成1,2,1,211nidMMMiiiii，。三种边界条件的三弯矩方程：（1）第一种边界条件：已知bfaf'',''。取)('',''0bfMafMn，这时方程组减少了两个未知量，变成只含n-1个未知量1,2,1niMi，的n-1个方程的方程组2,3,2,222111211102niMdMMdMMMMdMMnnnnnniiiiiiii，，用矩阵表示为nnnnnnnnnMdddMdMMMM112201112211222212222可用追赶法求解出)1,2,1(niMi，后，即得三条样条插值函数。（2）第二种边界条件，已知)(')('bfaf，，记)(''),(''0bfyafyn，则有')(',')('00nnyxSyxS，得'36'63-1101011101nnnnnnnnyhyyMhMhyhyyMhMh，，即nnndMMdMM2,21010，其中],,[6],,,[612100nnnnxxxfdxxxfd，得到方程组1,2,1,222111010nidMMdMMMdMMnnniiiiii，，用矩阵表示为nnnnnnddddMMMM1101101111212212，该方程组的系数矩阵是严格三对角占优矩阵，可用追赶法求解。（3）第三种边界条件：周期型边界条件。已知xfy是以0xxabTn为周期的周期函数，则由周期性可知，11110110,,,hhMMMMyyyynnnnn，，这时将点nx看成是内节点，则有nnnnnndMMM112，也即nnnnndMMM1112，其中nnnnnnhyyhyyhhd11116，方程组第1个方程为：121112dMMMn，所以方程组为1,3,2,222111111211nidMMMdMMMdMMMnnnnniiiiiin，，用矩阵表示为nnnnnnnnddddMMMM1211211122112222，显然系数矩阵为严格对角占优矩阵，可用LU分解法求解。2.程序框图计算步长得出三弯矩方程组解矩阵方程计算三次样条曲线的系数用插值函数计算出各点的插值计算值3.源程序functionx=mchase(A,d)%追赶法n=length(d);u=zeros(n,1);u(1)=A(1,1);fork=2:nl(k)=A(k,k-1)/u(k-1);u(k)=A(k,k)-l(k)*A(k-1,k);endy=zeros(n,1);y(1)=d(1);fori=2:ny(i)=d(i)-l(i)*y(i-1);endx=zeros(n,1);x(n)=y(n)/u(n);fori=n-1:-1:1x(i)=(y(i)-A(i,i+1)*x(i+1))/u(i);endxendfunctionT=mspline1(x0,y0,f21,f22,xx)%三次样条插值函数第一种边界条件%x0、y0分别为节点的横坐标和纵坐标；%f21为左端点的二阶导数值,f22为右端点的二阶导数值；xx为由插值点组成的向量n=length(x0)-1;%计算小区间数fori=1:nh(i)=x0(i+1)-x0(i);endfori=1:n-1mu(i)=h(i)/(h(i)+h(i+1));lamda(i)=1-mu(i);d(i)=6*((y0(i+2)-y0(i+1))/h(i+1)-(y0(i+1)-y0(i))/h(i))/(h(i)+h(i+1));endA=zeros(n-1);fori=1:n-2A(i+1,i)=mu(i+1);%次下对角线A(i,i+1)=lamda(i);%次上对角线A(i,i)=2;%主对角线endA(n-1,n-1)=2;dd=zeros(n-1,1);%右端列向量fori=2:n-2dd(i)=d(i);enddd(1)=d(1)-mu(1)*f21;dd(n-1)=d(n-1)-lamda(n-1)*f22;M=mchase(A,dd);%追赶法求解M值hmulamdaAddM=[f21,M',f22]'t=sym('t');a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);fori=1:na(i)=M(i)./(6*h(i));b(i)=M(i+1)./(6*h(i));W1(i)=b(i)-a(i);W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i));c(i)=y0(i)./h(i)-h(i).*M(i)/6;e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6;W3(i)=3*b(i).*x0(i).^2-3*a(i).*x0(i+1).^2+e(i)-c(i);W4(i)=a(i).*x0(i+1).^3-b(i).*x0(i).^3+c(i).*x0(i+1)-e(i).*x0(i);Si(t)=W1(i).*t^3+W2(i).*t^2+W3(i).*t+W4(i)%每个小区间的三次样条差值函数表达式endm=length(xx);T=zeros(m,1);fork=1:mforj=1:nif((xx(k)=x0(j))&(xx(k)x0(j+1)))T(k)=W1(j).*(xx(k).^3)+W2(j).*(xx(k).^2)+W3(j).*xx(k)+W4(j);endendendTEndfunctionT=mspline2(x0,y0,f11,f12,xx)%三次样条插值函数第二种边界条件%x0、y0分别为节点的横坐标和纵坐标；%f11为左端点的二阶导数值,f12为右端点的二阶导数值；xx为由插值点组成的向量n=length(x0)-1;%计算小区间数fori=1:nh(i)=x0(i+1)-x0(i);endfori=1:n-1mu(i)=h(i)/(h(i)+h(i+1));lamda(i)=1-mu(i);d(i)=6*((y0(i+2)-y0(i+1))/h(i+1)-(y0(i+1)-y0(i))/h(i))/(h(i)+h(i+1));endA=zeros(n+1);%系数矩阵dd=zeros(n+1,1);%右端列向量fork=2:nA(k,k)=2;%主对角线元素A(k,k-1)=mu(k-1);%次下对角线元素A(k,k+1)=lamda(k-1);%次上对角线元素endA(1,1)=2;A(1,2)=1;A(n+1,n+1)=2;A(n+1,n)=1;dd(1)=6*((y0(2)-y0(1))/h(1)-f11)/h(1);dd(n+1)=6*(f12-(y0(n+1)-y0(n))/h(n))/h(n);fork=2:ndd(k)=d(k-1);endM=mchase(A,dd);%追赶法求解M值hmulamdaAddMt=sym('t');a=zeros(n,1);b=zeros(n,1);c=zeros(n,1);e=zeros(n,1);fori=1:na(i)=M(i)./(6*h(i));b(i)=M(i+1)./(6*h(i));W1(i)=b(i)-a(i);W2(i)=3*(a(i).*x0(i+1)-b(i).*x0(i));c(i)=y0(i)./h(i)-h(i).*M(i)/6;e(i)=y0(i+1)./h(i)-h(i).*M(i+1)/6;W3(i)=3*b(i).*x0(i).^2-3*a(i).*x0(i+1).^2+e(i)-c(i);W4(i)=a(i).*x0(i+1).^3-b(i).*x0(i).^3+c(i).*x0(i+1)-e(i).*x0(i);Si(t)=W1(i).*t^3+W2(i).*t^2+W3(i).*t+W4(i)%每个小区间的三次样条差值函数表达式endm=length(xx);T=zeros(m,1);fork=1:mforj=1:nif((xx(k)=x0(j))&(xx(k)x0(j+1)))T(k)=W1(j).*(xx(k).^3)+W2(j).*(xx(k).^2)+W3(j).*xx(k)+W4(j);endendendTend%计算实习，第一种边界条件clear;x=-1:0.2:1;%输入节点横坐标y=ones(1,11)./(ones(1,11)+25*x.^2)%计算节点纵坐标t=sym('t');f=1/(1+25*t^2);%