用三次樣條插值求離散點(diǎn)斜率 matlab程序

shouce

用三次樣條插值求斜率 三次樣條插值的matlabt程序 function x=followup(a,b,c,d)n=length(d); 0 Z, r% j2 O- u5 H- j6 Ea(1)=0; 9 R  D- c* G8 W. m0 N6 [%“追”的過程 ; n" i  z( k# B8 O" L: `, I, P4 WL(1)=b(1); / D/ c% }; W& _y(1)=d(1)/L(1); 2 i- V# T0 ], m2 a' p, q7 Zu(1)=c(1)/L(1); / Z% E+ r) R: J/ D+ _, [/ C& Xfor i=2n-1)     L(i)=b(i)-a(i)*u(i-1);     y(i)=(d(i)-y(i-1)*a(i))/L(i);     u(i)=c(i)/L(i);   ?+ a1 d0 v  f5 k6 i) [8 jend 7 B" y! g6 c. e+ KL(n)=b(n)-a(n)*u(n-1); - r$ b' {; f2 v+ [y(n)=(d(n)-y(n-1)*a(n))/L(n); . O7 d7 X/ Q& c/ {  u0 X%“趕”的過程 x(n)=y(n); / z, s9 I8 p' _3 ?* A- Xfor i=(n-1):-1:1 * v9 P" D0 \4 e6 I/ j    x(i)=y(i)-u(i)*x(i+1); end 7 N3 D: j; j- s, g function[s,y0]=spline3 (x,y,x0) %x,y為數(shù)表x0為插值點(diǎn)s表示插值函數(shù)y0為x0對(duì)應(yīng)的插值函數(shù)值 : \1 H' }3 A# Fsyms t n=length(x); %得出n ! s, h8 H# z/ a' }2 ^% H* Ofor i=1:n-1; 1 M8 {4 e+ ]8 {  t! p, \) \    h(i)=x(i+1)-x(i); ! j8 Y3 V; A- `8 F2 n8 Bend , k- n# p: u5 m$ P3 hfor i=2:n-1; # A+ u% `# @/ l6 k- x+ m6 ]    lamda(i)=h(i)/(h(i-1)+h(i));     miu(i)=1-lamda(i);   B0 o; f- }6 B6 I9 U# S    g(i)=3*(lamda(i)*((y(i)-y(i-1))/h(i-1))+miu(i)*((y(i+1)-y(i))/h(i))); end g(1)=3*((y(2)-y(1))/h(1)); g(n)=3*((y(n)-y(n-1))/h(n-1)); $ b" y( K/ f  c- R! `6 }+ w0 l%前邊求出lamda miu和g從而可以確定系數(shù)矩陣 : |. Z5 L8 v8 J8 f# R/ J& e2 mmiu(1)=1;   U9 V  U, J1 X' q9 V2 Bmiu(4)=0; . Z- X9 f/ b6 H* Y& e! Olamda(n)=1; lamda(1)=0; - ?0 q* n6 m5 z$ }+ x%根據(jù)第二邊界條件補(bǔ)充兩個(gè)lamda和miu的值 for i=1:n - I, A7 S$ }' l    beta(i)=2; end 4 y" c- d& U" E8 pm=followup(lamda,beta,miu,g) %解出m的值從而可確定st st為各段的插值多項(xiàng)式 8 O3 K7 ?% N7 R) p) ]# \/ ]for i=1:n-1 ( G: q8 c8 s  M" X    st(i)=(t-x(i+1))^2*(h(i)+2*(t-x(i)))*y(i)/(h(i)^3)... & s1 ~8 Q/ H8 p( Z    +(t-x(i))^2*(h(i)+2*(x(i+1)-t))*y(i+1)/(h(i)^3)...   G# c6 w% _8 g* X    +(t-x(i))^2*(t-x(i+1))*m(i+1)/(h(i)^2)...     +(t-x(i+1))^2*(t-x(i))*m(i)/(h(i)^2); ! H0 a' E+ u/ x: i2 uend %得到插值的結(jié)果各段的t的表達(dá)式 %接下來要將插值點(diǎn)x0代入首先確定x0所在的插值區(qū)間 for i=1:n-1 , y  K7 s$ g" O9 l8 }    if (x(i)>x0) & d+ Y% f' b" Z1 t& j        in=i;     end % T3 v6 b+ S, H* I$ M) p  Eend 1 X7 G4 O! ?4 [, a8 hs=st(in); , n  y  e: ?, X0 q  j0 Xs=expand(s); s=collect(s,'t'); + p. o6 m, t9 l9 Z8 w2 v2 xy0=subs(s,'t',x0) %s是插值多項(xiàng)式y(tǒng)0是插值點(diǎn)的函數(shù)值 在matlab中輸入 x=[1 2 4 5]； y=[1 3 4 2]； spline3(x,y,2) 9 w. a( c3 I7 m會(huì)得到各點(diǎn)的斜率 + f7 L: X+ _: t3 H