资源描述:
《数值分析教教案7》由会员上传分享,免费在线阅读,更多相关内容在工程资料-天天文库。
1、1.7.4三次样条插值的MATLAB实现实例【例「21】在区间的两个端点处给出条件,称为边界条件。给定两端点处的导数值s’S)=%,s'(b)=元,然后根据下面的数据点求出其三次样条插值多项式,并计算当兀=3・5时y的值。Xi245678y0.840.910.14-0.76-0.96-0.280.660.99/y0.540341455解:在MATLAB编辑窗口中输入以下命令:x=l:&y二[0.840.910.14-0.76-0.96-0.280.66,0.99][f,fO]=ThrSamplel(x,y,0.54,-0.15,3-5)程序运行后得到如下结果。f=(-3
2、8/25*t+133/25)*(t-5)A2+(-264/25+48/25*t)*(t-4)A2+(-944刀14555逍+37788/14555)*(5-t)A2-(25502967261225075/18014398509481984-5100593452245015/18014398509481984*t)*(t-4)A2f0二-0.9164程序运行时调用的三次样条函数如下(给定两端点处的一阶导数值)。function[f,fO]二ThrSamplel(x,y,y_l,y_N,xO)symst;f=0・0;fO二0.0;if(length(x)==length(y
3、))n=length(x);elsedisp('x和y的维数不相等!’);return;%维数检查endfori=l:nif(x(i)<=x0)&&(x(i+l)>=xO)index=i;break;endend%找到xO所在区间A=diag(2*ones(l,n));%求解m的系数矩阵u=zeros(n-2,l);lamda=zeros(n-l,l);c二zeros(n,l);fori=2:n-lu(i-l)=(x(i)-x(i-l))/(x(i+l)-x(i-l));lamda(i)二(x(i+l)-x(i))/(x(i+l)-x(i-l));c(i)=3*lamd
4、a(i)*(y⑴-y(i-1))/(x(i)-x(i-1))+...3*u(i・1)*(y(i+l)-y(i))/(x(i+l)-x(i));endMl)=2*y_l;c(n)=2*y_N;m=followup(A,c);%用追赶法求解方程组h=x(index+l)-x(index);%xO所在区间长度f=y(index)*(2*(t-x(index))+h)*(t-x(index+l))八2/h/h/h+...y(index+l)*(2*(x(index+l)-t)+h)*(t-x(index))A2/h/h/h+...m(index)*(t-x(index))*(x
5、(index+l)-t)A2/h/h-…m(index+l)*(x(index+l)-t)*(t-x(index))A2/h/h;%xO所在区间的插值函数fO=subs(f,f,x0);%x0处的插值下面函数是ThrSamplel(x,y,y_l,y_N,xO)函数调用运行时调用的函数。functionx=followup(A,b)n=rank(A);for(i=l:n)if(A(i,i)==0)disp(!Error:对角有元素为0!J;return;endend;d=ones(n,l);a=ones(n-l,l);c=ones(n-l);for(i=l:n-l)a(
6、i,l)=A(i+l,i);c(i,l)=A(i,i+l);d(i,l)=A(i,i);endd(n,l)=A(n,n);for(i=2:n)d(i,l)=d(i,l)-(a(i-l,l)/d(i-l,l))*c(i-l,l);b(i,l)=b(i,l)-(a(i-l,l)/d(i-l,l)rb(i-l,l);endx(n,l)=b(n,l)/d(n,l);for(i=(n-l):-l:l)x(i,l)=(b(i,l)-c(i,l)*x(i+l,l))/d(i,l);end【例1-22]在区间的两个端点处给出条件,称为边界条件。给定两端处的二阶导数值S(/)二二£,然后
7、根据下面的数据点求出其三次样条插值多项式,并计算当%=3・5时y的值。X12345678y0.840.910.14-0.76-0.96-0.280.660.99y-0.8415-0.9894解:在MATLAB编辑窗口中输入以下命令:x=l:8;y=[0.840.910.14-0.76-0.96-0.280.660.99];[f,f0]=ThrSample2(x,y,-0.84,-0.99,3.5)程序运行后得到如下结果。(7/25*t・7/10)*(t・4)八2+(-171/25+38/25*t)*(t・3)八2+(-4455276049