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

这一题是得到数据点(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的范围。

    最后要输出表达式的话这个还是有点复杂的：使用以下函数。9 p) H: _( x7 v; @4 V8 Y' N
       pp=interp1(ax,ay,'spline','pp');; x* N9 k0 F0 w2 u% Y0 }; W4 u
       breaks=pp.breaks; %breaks的第i和i+1个元素为m和n
- y9 e. I. K. v& o       coefs=pp.coefs; %假设coefs的第i行为a b c d，
# \$ ~" Q" k! R+ d& ^1 o    接着再用一个循环得出每个表达式输出各个表达式即可。

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

   % use natural cubic spline to interpolate data point
. ^  u! i$ w9 T6 m; \& l& {! j' t   % a、(0,3),(1,5),(2,4),(3,1)
% y& u) @8 @. S8 c3 Z! Q3 C   % b、(-1,3),(0,5),(3,1),(4,1),(5,1)
' \( E6 V1 \/ h( E# sfunction page_178_1_natural_spline_script
/ A3 p8 ~& D9 ?& l" E% L- {; sax = [0 1 2 3];
6 Q! ]& A  v3 bay = [3 5 4 1];  %对数据点(0,3),(1,5),(2,4),(3,1)进行三次样条建模，并输出表达式和图像
& z# A. D" |2 e5 x$ `" Nbx = [-1 0 3 4 5];6 v5 u5 v8 t- e& o
by = [3 5 1 1 1];%对数据点(-1,3),(0,5),(3,1),(4,1),(5,1)进行三次样条建模，并输出表达式和图像
" I- }9 ]/ y2 r0 f) }3 `csa = spline(ax,[0 ay 0]);
$ P4 G: ~  A; z) k# R! m+ G( c/ qxxa = linspace(0,3,1000);
6 U1 [  I, C) w$ Csubplot(1,2,1);
. i) k! X. B) i$ hplot(ax,ay,'o',xxa,ppval(csa,xxa),'-');
' e) V2 v+ j: j! e! L. ~( j' vxlabel('a x 0~3');
4 u* K, y  P5 d( D: y0 X6 ]ylabel('a y');! O+ i3 L1 T6 p; C
title('equation a');
: S: z  b3 F7 ~; a& hcsb = spline(bx,[0 by 0]);
1 `, _/ {2 i* A- e0 F( Z3 kxxb = linspace(-1,5,10000);
2 W# c8 S# z! V; \8 O5 Z  G  M# ksubplot(1,2,2);) R6 N2 [* G0 O( o  q3 b
plot(bx,by,'o',xxb,ppval(csb,xxb),'-');: P* C7 f. y/ G/ I5 {
xlabel('b x -1~5');
7 S7 e- k. h, z$ v5 H% W# tylabel('b y');4 ?( e2 |8 ^, K' F5 n
title('equation b');

pp=interp1(ax,ay,'spline','pp');
( b3 G& B0 v& H7 w/ Mbreaks=pp.breaks; %breaks的第i和i+1个元素为m和n
- @6 u5 H9 P. q: m) D  @3 Pcoefs=pp.coefs; %假设coefs的第i行为a b c d，
2 x5 q. x) ~2 C) a0 ]     %breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是
$ Y# b8 A# [% s  ~     %a(x-m)^3+b(x-m)^2+c(x-m)+d
+ {& b4 o& q, ]syms x
0 K6 e$ T8 J% o2 o8 u5 M+ n0 `. xdisp('对于数据点(0,3),(1,5),(2,4),(3,1)的三次样条表达式为：');) R$ d; v. a1 S- ]
for i = 1:3
6 P' N( q" w4 z2 G- r, b. a1 `    y = coefs(i,1)*((x - breaks(i))^3) + coefs(i,2)*((x - breaks(i))^2) + coefs(i,3)*((x - breaks

    (i))) + coefs(i,4), ?5 k" e: V2 S$ D' M9 D
end
" A, x! s8 y3 [2 t3 o( ^- N1 Bppb=interp1(bx,by,'spline','pp');
, w' c+ n: h' C, xbreaksb=ppb.breaks; %breaks的第i和i+1个元素为m和n7 a4 C9 J- m: D0 r( {
coefsb=ppb.coefs; %假设coefs的第i行为a b c d，
7 E5 Z2 n# T' }" z3 l4 I5 d%breaks的第i和i+1个元素为m和n那么在区间[m,n]的函数表达式就是
. _% C3 O* ^& |, n%a(x-m)^3+b(x-m)^2+c(x-m)+d
+ F" Y( E- H/ j1 d) Lsyms x9 W) F+ m1 z- A" P+ [' @
disp('对于数据点(-1,3),(0,5),(3,1),(4,1),(5,1)的三次样条表达式为：');
" i& z6 g/ m$ E, v5 A; x+ xfor i = 1:4
* _3 a! p* b% }% 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)! G* G; i( i# w7 h0 Y9 \* |* h- N
end