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

shouce

1 s2 ~1 k6 [* @, V$ P" ]( L
! h% |# f2 r# j, |

用三次样条插值求斜率 三次样条插值的matlabt程序 function x=followup(a,b,c,d)n=length(d); a(1)=0; %“追”的过程 L(1)=b(1); y(1)=d(1)/L(1); ' |9 z8 ~: a$ H* g. Q' Ou(1)=c(1)/L(1); * O' A+ I# `8 M+ y  S. F/ P# T" t7 Hfor i=2n-1) & h6 K% g  j2 q: j. c    L(i)=b(i)-a(i)*u(i-1);     y(i)=(d(i)-y(i-1)*a(i))/L(i); 5 {: F( [0 L, V4 h+ `0 q( r    u(i)=c(i)/L(i); : `) \: P; N1 |) w4 l% Pend L(n)=b(n)-a(n)*u(n-1);   z6 _% w% P  H( ?/ v/ f- J$ Qy(n)=(d(n)-y(n-1)*a(n))/L(n); %“赶”的过程 6 G# p3 u9 E9 e; f! u9 Kx(n)=y(n); 7 }" L$ `% D5 ^for i=(n-1):-1:1 ' y& y- |+ p' ~. N: L  U    x(i)=y(i)-u(i)*x(i+1); 9 |# s: _- x" Y5 h4 s. N- Pend 9 l4 D8 B" g+ P/ w- O' Q 8 ]  `4 L" k7 _$ U function[s,y0]=spline3 (x,y,x0) %x,y为数表x0为插值点s表示插值函数y0为x0对应的插值函数值 syms t / d0 E' W& K' G, Wn=length(x); %得出n % O' D5 h; R. x( Xfor i=1:n-1; 4 o2 t& @6 q5 i/ X, j( O7 M1 n    h(i)=x(i+1)-x(i); end . F6 H- B& f" `3 q: t9 Z- Nfor i=2:n-1; ) `: ?6 F9 z/ b1 J    lamda(i)=h(i)/(h(i-1)+h(i));     miu(i)=1-lamda(i);     g(i)=3*(lamda(i)*((y(i)-y(i-1))/h(i-1))+miu(i)*((y(i+1)-y(i))/h(i))); " G0 D, R4 s6 ^+ n+ ~9 x! l* K  rend 7 E/ M* t2 P/ W% N" J. g4 ~3 [g(1)=3*((y(2)-y(1))/h(1)); 7 P, B+ S. G5 t$ V4 W) S5 E* H) ?* hg(n)=3*((y(n)-y(n-1))/h(n-1)); %前边求出lamda miu和g从而可以确定系数矩阵 + Z% G+ h3 d, ^! l: \5 Cmiu(1)=1; miu(4)=0; ' V& Y3 C' ^( y& k% T1 Llamda(n)=1; 3 l) n8 U0 G. Y$ s4 F& Q8 l5 }" ilamda(1)=0; " p8 b/ K/ l1 o) _%根据第二边界条件补充两个lamda和miu的值 for i=1:n : O6 S$ x* T3 ^8 A, p8 f3 y    beta(i)=2; end 8 J" Z2 x% i" cm=followup(lamda,beta,miu,g) %解出m的值从而可确定st st为各段的插值多项式 for 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); end   `. b+ I2 X3 z" }3 }%得到插值的结果各段的t的表达式 %接下来要将插值点x0代入首先确定x0所在的插值区间 for i=1:n-1     if (x(i)>x0) $ v/ e3 r* S, e' ?' A        in=i; . [; c8 r8 W7 I' M1 g) D    end + L+ }* b( n4 P, Nend . [1 `! v9 s. r! H. As=st(in); & c) @5 {+ M) ns=expand(s); 4 l, J, V/ O3 y8 K* ks=collect(s,'t'); y0=subs(s,'t',x0) ; G/ ^$ H) L' f%s是插值多项式y0是插值点的函数值 2 R" I0 r' \9 {7 `* Q % T: X. N% h2 I8 Y/ k6 j; Z7 J 在matlab中输入 x=[1 2 4 5]； y=[1 3 4 2]； spline3(x,y,2) 7 U5 B3 v7 d5 n) ~1 f8 N会得到各点的斜率 " u4 X$ \! b8 ]9 j! |4 D0 N! ~

. X% ^) Z% A, x7 h