这一题是得到数据点(0,3),(1,5),(2,4),(3,1)并得到它的三次样条表达式和画出三次样条后的图图形。
以及对数据点(-1,3),(0,5),(3,1),(4,1),(5,1)并得到它的三次样条表达式和画出三次样条后的图图形。
用函数spline可以直接得到,都是如果是要求自然三次样条呢?那就可以在数组y的左右两侧添0。如:
csa = spline(ax,[0 ay 0]);再用xxa = linspace(0,3,1000);plot(xxa,ppval(csa,xxa));进行绘图。linespace是将0~3分成1000份,然后ppval是求三次样条在不同的xxa上的值。MATLAB中的插值函数为interp1,其调用格式为: yi= interp1(x,y,xi,'method'),其中x,y为插值点,yi为在被插值点xi处的插值结果;x,y为向量, 'method'表示采用的插值方法,MATLAB提供的插值方法有几种: 'method'是最邻近插值, 'linear'线性插值; 'spline'三次样条插值; 'cubic'立方插值.缺省时表示线性插值。注意:所有的插值方法都要求x是单调的,并且xi不能够超过x的范围。
最后要输出表达式的话这个还是有点复杂的:使用以下函数。 pp=interp1(ax,ay,'spline','pp'); breaks=pp.breaks; %breaks的第i和i+1个元素为m和n coefs=pp.coefs; %假设coefs的第i行为a b c d, 接着再用一个循环得出每个表达式输出各个表达式即可。
一下是我的代码,写的有一点粗糙,希望别见怪啊!
% use natural cubic spline to interpolate data point % a、(0,3),(1,5),(2,4),(3,1) % b、(-1,3),(0,5),(3,1),(4,1),(5,1) function page_178_1_natural_spline_script ax = [0 1 2 3]; ay = [3 5 4 1]; %对数据点(0,3),(1,5),(2,4),(3,1)进行三次样条建模,并输出表达式和图像 bx = [-1 0 3 4 5]; by = [3 5 1 1 1];%对数据点(-1,3),(0,5),(3,1),(4,1),(5,1)进行三次样条建模,并输出表达式和图像 csa = spline(ax,[0 ay 0]); xxa = linspace(0,3,1000); subplot(1,2,1); plot(ax,ay,'o',xxa,ppval(csa,xxa),'-'); xlabel('a x 0~3'); ylabel('a y'); title('equation a'); csb = spline(bx,[0 by 0]); xxb = linspace(-1,5,10000); subplot(1,2,2); plot(bx,by,'o',xxb,ppval(csb,xxb),'-'); xlabel('b x -1~5'); ylabel('b y'); title('equation b');
pp=interp1(ax,ay,'spline','pp'); breaks=pp.breaks; %breaks的第i和i+1个元素为m和n coefs=pp.coefs; %假设coefs的第i行为a b c d, %breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是 %a(x-m)^3+b(x-m)^2+c(x-m)+d syms x disp('对于数据点(0,3),(1,5),(2,4),(3,1)的三次样条表达式为:'); for i = 1:3 y = coefs(i,1)*((x - breaks(i))^3) + coefs(i,2)*((x - breaks(i))^2) + coefs(i,3)*((x - breaks
(i))) + coefs(i,4) end ppb=interp1(bx,by,'spline','pp'); breaksb=ppb.breaks; %breaks的第i和i+1个元素为m和n coefsb=ppb.coefs; %假设coefs的第i行为a b c d, %breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是 %a(x-m)^3+b(x-m)^2+c(x-m)+d syms x disp('对于数据点(-1,3),(0,5),(3,1),(4,1),(5,1)的三次样条表达式为:'); for i = 1:4 y = coefsb(i,1)*((x - breaksb(i))^3) + coefsb(i,2)*((x - breaksb(i))^2) + coefsb(i,3)*((x -
breaksb(i))) + coefsb(i,4) end