1. 程式人生 > >Matlab之三次樣條畫圖和表示式

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的範圍。

    最後要輸出表達式的話這個還是有點複雜的:使用以下函式。
       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