Matlab之三次样条画图和表达式

haidaowang

这一题是得到数据点(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
1 ?( N* T5 Y& Q6 R9 w       coefs=pp.coefs; %假设coefs的第i行为a b c d，
    接着再用一个循环得出每个表达式输出各个表达式即可。

    一下是我的代码，写的有一点粗糙，希望别见怪啊！

   % use natural cubic spline to interpolate data point
9 E+ _  U2 `: x2 K% b4 w3 u   % a、(0,3),(1,5),(2,4),(3,1)
$ `1 Y0 C- B  y  B   % b、(-1,3),(0,5),(3,1),(4,1),(5,1)
& x( y# Q3 J7 P; I5 r- n& [6 ofunction page_178_1_natural_spline_script
ax = [0 1 2 3];
ay = [3 5 4 1];  %对数据点(0,3),(1,5),(2,4),(3,1)进行三次样条建模，并输出表达式和图像
+ f3 @+ F/ F$ {  G# ]bx = [-1 0 3 4 5];
2 Q1 ^1 Q; y* |3 F# D3 Iby = [3 5 1 1 1];%对数据点(-1,3),(0,5),(3,1),(4,1),(5,1)进行三次样条建模，并输出表达式和图像
csa = spline(ax,[0 ay 0]);
% }8 a5 e  c6 `( _& Zxxa = linspace(0,3,1000);
! P/ W0 K' g% T% p6 lsubplot(1,2,1);
+ ?! c: c- A* B- R7 I4 O& S( c8 dplot(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);
' i/ C! I; u/ w9 Rplot(bx,by,'o',xxb,ppval(csb,xxb),'-');
4 m. n- o4 H6 G; K& X) `/ rxlabel('b x -1~5');
/ t" ?* \/ b  D" p( ]ylabel('b y');
1 n! ^/ O6 p  O( u: ftitle('equation b');

pp=interp1(ax,ay,'spline','pp');
8 P1 r) Y- ^. {breaks=pp.breaks; %breaks的第i和i+1个元素为m和n
- k; g5 [) x6 U8 k  |# r) fcoefs=pp.coefs; %假设coefs的第i行为a b c d，
' R) V* A& u  X5 ?     %breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是
- d0 F( J5 U/ C$ Z) ^     %a(x-m)^3+b(x-m)^2+c(x-m)+d
syms x
; C# P, m5 W6 s* f( ~+ _1 vdisp('对于数据点(0,3),(1,5),(2,4),(3,1)的三次样条表达式为：');
* Q2 P  S; g- `# t5 ifor i = 1:3
3 m3 f1 U( \6 A* ]" ^2 ?    y = coefs(i,1)*((x - breaks(i))^3) + coefs(i,2)*((x - breaks(i))^2) + coefs(i,3)*((x - breaks

    (i))) + coefs(i,4)
' r+ P: a% V6 Z& e8 K: X  c2 Tend
( f. s, ~2 H1 f: V' `4 @$ Eppb=interp1(bx,by,'spline','pp');
breaksb=ppb.breaks; %breaks的第i和i+1个元素为m和n
! j' {% e( V3 m' ^) V! icoefsb=ppb.coefs; %假设coefs的第i行为a b c d，
( f! \3 J8 g6 }, {%breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是
0 V* w1 X/ J6 B2 J7 c%a(x-m)^3+b(x-m)^2+c(x-m)+d
5 _* f% ^6 O- s! e2 s8 u3 v% psyms x
; z" c" N0 g1 l$ [disp('对于数据点(-1,3),(0,5),(3,1),(4,1),(5,1)的三次样条表达式为：');
9 d/ s: q7 k4 c& A4 S7 hfor i = 1:4
7 q2 [  F# x6 W8 j; F$ {( O6 F/ d% [' U   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

' N6 ]' O( |& _3 o. o5 n& J% d5 s

BeardRen

学习