用三次样条插值求离散点斜率 matlab程序

shouce

( ?' E* L% P+ w- S; F

用三次样条插值求斜率 三次样条插值的matlabt程序 ) l- h) g" g7 ^* ?  Y function x=followup(a,b,c,d)n=length(d); 8 I. k/ }* \& I) |1 na(1)=0; %“追”的过程 L(1)=b(1); - _2 |$ m* i  w2 a3 U% K  \y(1)=d(1)/L(1); u(1)=c(1)/L(1); ' |* O+ @1 g( t" [for i=2n-1) 3 ^4 |  h3 F" s! m4 P- a* D    L(i)=b(i)-a(i)*u(i-1);     y(i)=(d(i)-y(i-1)*a(i))/L(i); + D: |2 M8 |" I3 z2 G    u(i)=c(i)/L(i); end L(n)=b(n)-a(n)*u(n-1); & j7 z3 @$ l2 p; l* ~# `# yy(n)=(d(n)-y(n-1)*a(n))/L(n); %“赶”的过程 / k( S0 o$ m" }' jx(n)=y(n); for i=(n-1):-1:1 ) }9 u6 Y: F0 K  v8 X1 D; f! ^    x(i)=y(i)-u(i)*x(i+1); % X% g. x. Y1 u! n4 A& s5 E$ lend 0 x$ h$ _5 Y/ L" i- c% L ) o% Z3 E! q% g' s( P! L. }4 R & W# v! a) ^3 o- ` function[s,y0]=spline3 (x,y,x0) 9 M# h: _; _8 i, J%x,y为数表x0为插值点s表示插值函数y0为x0对应的插值函数值 syms t ; f* ?6 C% f0 e* t9 C- [& M% kn=length(x); %得出n for i=1:n-1;     h(i)=x(i+1)-x(i); end for i=2:n-1;     lamda(i)=h(i)/(h(i-1)+h(i));     miu(i)=1-lamda(i); $ L% S# }) u) f- C' v$ \# i' z    g(i)=3*(lamda(i)*((y(i)-y(i-1))/h(i-1))+miu(i)*((y(i+1)-y(i))/h(i))); end 7 b0 X! F5 I% n3 E! J& T$ t# w! |g(1)=3*((y(2)-y(1))/h(1)); g(n)=3*((y(n)-y(n-1))/h(n-1)); 9 R4 e( q& |+ I4 D8 i9 ]* I1 Y%前边求出lamda miu和g从而可以确定系数矩阵 % _" b0 T" }% f9 Z2 m- ]( Y! X7 E. W* Gmiu(1)=1; % }& Z- d& S9 P0 Lmiu(4)=0; lamda(n)=1; lamda(1)=0; 0 C: B6 z; A1 ]8 f3 K; d/ d%根据第二边界条件补充两个lamda和miu的值 for i=1:n 5 Y( E  X& ]  M* N    beta(i)=2; 6 X6 a" [+ x5 d+ `: N2 w! }end m=followup(lamda,beta,miu,g) %解出m的值从而可确定st st为各段的插值多项式 4 H' X  A4 D% C) Cfor i=1:n-1     st(i)=(t-x(i+1))^2*(h(i)+2*(t-x(i)))*y(i)/(h(i)^3)...     +(t-x(i))^2*(h(i)+2*(x(i+1)-t))*y(i+1)/(h(i)^3)...     +(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); 7 S. N; b# F( E& C) Kend %得到插值的结果各段的t的表达式 %接下来要将插值点x0代入首先确定x0所在的插值区间   k/ ^( W7 q0 I' E) lfor i=1:n-1     if (x(i)>x0)         in=i;     end end 5 B+ O! J) L! b0 `8 D8 m$ Us=st(in); s=expand(s); 6 k( T2 p4 r) O4 D$ |0 s3 u# Hs=collect(s,'t'); 4 H; ~$ Y( x9 A- A9 Ey0=subs(s,'t',x0) %s是插值多项式y0是插值点的函数值 3 y" o+ X" Z* x( n; I" S3 D % ^! @! x1 a  z7 ?( ]7 n" v9 l 在matlab中输入 x=[1 2 4 5]； y=[1 3 4 2]； spline3(x,y,2) 会得到各点的斜率 / \- |, c5 X4 ` : t2 o: l: E8 ?