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

三次样条插值的matlabt程序

: _0 {& i) ^* K, b- i: tfunction x=followup(a,b,c,d)n=length(d);* i0 s$ G& t6 T$ ?8 a7 m

a(1)=0;7 Q8 c5 R. `/ N# b

%“追”的过程( d$ ^7 g7 F" W% j) |( T2 ?0 Y) e; X

L(1)=b(1);7 B% ~( G/ X+ G: U3 G5 z4 }

y(1)=d(1)/L(1);

$ |# `/ f+ {) I3 eu(1)=c(1)/L(1);8 ?8 I! W8 V2 S  h% f' |. n7 I) j6 c

for i=2

n-1)

0 H; K  ]6 O, y$ ]$ M' A$ W6 rL(i)=b(i)-a(i)*u(i-1);& Z# Z8 `  B. L; B

y(i)=(d(i)-y(i-1)*a(i))/L(i);5 d' Y) f. E/ F+ c  `5 H; _; l

u(i)=c(i)/L(i);

5 ^& q: w$ f0 r! T% q. _, _0 H1 uend- Y0 g+ n( M" |9 c

L(n)=b(n)-a(n)*u(n-1);4 G% A+ ^- O$ n: z9 O

y(n)=(d(n)-y(n-1)*a(n))/L(n);/ B5 l* p: p$ Y% X

%“赶”的过程0 }$ C5 u- _( b0 ], t

x(n)=y(n);" a" b: y: l* B3 X) k& B9 [

for i=(n-1):-1:1

5 O1 K/ }9 l0 h6 @( ?/ A  P$ mx(i)=y(i)-u(i)*x(i+1);; T3 h# w. m8 w# L, ?5 _! o+ H! B

end

, Z1 \5 c, z6 z" J6 g* k- ?' |" `1 M1 t/ a; i. Z3 q

& X; N- ^! F) p1 f% X- {; mfunction[s,y0]=spline3 (x,y,x0)" [% k6 ?! [% Z" K4 ~

%x,y为数表x0为插值点s表示插值函数y0为x0对应的插值函数值}& \* z: U' B4 B6 W- ^8 Q( M

syms t

* }/ @8 C3 H1 X  g6 j; a' z9 En=length(x);

0 L  ]) n7 N6 w3 H- t  M" w%得出n

0 \* L1 f, a" h1 V( Ufor i=1:n-1;

' g2 K# G; B' }, t( ?h(i)=x(i+1)-x(i);

4 i( p; P5 _0 D+ bend& o( n0 o  j3 E' c' X

for i=2:n-1;

8 J' M) J8 [* _+ O! Wlamda(i)=h(i)/(h(i-1)+h(i));

0 ~  o9 [, o1 m1 Dmiu(i)=1-lamda(i);

3 @- k5 _3 M! e8 [  k% ~0 F' ~* rg(i)=3*(lamda(i)*((y(i)-y(i-1))/h(i-1))+miu(i)*((y(i+1)-y(i))/h(i)));

. ?$ b! r0 |3 t( rend

% p5 y" a& f  l: p. |  J" H1 og(1)=3*((y(2)-y(1))/h(1));) |; y9 v8 M! E1 a5 `) t

g(n)=3*((y(n)-y(n-1))/h(n-1));

4 G! W4 C2 P/ z%前边求出lamda miu和g从而可以确定系数矩阵5 H6 T, D3 N+ _3 R( |( a( {

miu(1)=1;4 H1 f" q5 A; N3 `- \

miu(4)=0;# \& [9 E8 `8 Z

lamda(n)=1;. O/ M, V+ R8 E) X' k

lamda(1)=0;

! ?/ L0 U4 P% Z+ q8 o  G$ K, B/ r%根据第二边界条件补充两个lamda和miu的值

# i! p8 V! E& B/ Dfor i=1:n

r) g1 H( Z' f# N/ hbeta(i)=2;" Q3 t2 S& p& y7 s3 L& v! M4 f

end

) @( ]/ A% ?5 \" X7 Sm=followup(lamda,beta,miu,g)r! A8 c) ^( N6 ?

%解出m的值从而可确定st st为各段的插值多项式

( p# H. z/ \) O) Bfor i=1:n-1

?; w! j! e1 g3 G' b) qst(i)=(t-x(i+1))^2*(h(i)+2*(t-x(i)))*y(i)/(h(i)^3)...

* t+ N0 C4 `0 L' [+(t-x(i))^2*(h(i)+2*(x(i+1)-t))*y(i+1)/(h(i)^3)...

3 |) A/ x# E" W9 q+ z6 ?- u7 k+(t-x(i))^2*(t-x(i+1))*m(i+1)/(h(i)^2)...

! T, c5 R5 J# M( z+(t-x(i+1))^2*(t-x(i))*m(i)/(h(i)^2);& l# c7 t$ Y8 U- R

end

" s* q2 b) ?5 s1 z& B* c%得到插值的结果各段的t的表达式! N# h2 }7 ^- U& D2 F: D

%接下来要将插值点x0代入首先确定x0所在的插值区间% X+ ]0 ?: J) f8 m9 J! }

for i=1:n-1

; P2 V3 t5 K: S; \0 T, T7 Rif (x(i)>x0)1 u9 t/ r- K) M# z

in=i;# `1 w& f9 [; F4 j- ?2 E1 o

end

/ c( \3 E, T& h- pend

) W! {( W, u. R8 Ts=st(in);

& u$ v7 J9 [6 ms=expand(s);1 @2 Z% w1 {$ U( w/ V: V

s=collect(s,'t');

, U9 z8 N" X) j8 J" ty0=subs(s,'t',x0)5 H0 r9 f) a+ `

%s是插值多项式y0是插值点的函数值+ @: [( v2 P' O1 Z% r. p. n4 A

4 s* Y6 E, \8 B& d2 e/ d' e

$ |6 ]/ u: ]) M+ u5 B8 i. c" k在matlab中输入

' |6 g: a: K0 M- [

x=[1 2 4 5]；

y=[1 3 4 2]；

spline3(x,y,2)

+ \( K8 e9 ]' L' m会得到各点的斜率/ m* [7 h  f0 s8 `% I1 A

! m. A# j5 T. V/ [! n8 a6 I" ^' s( \, k2 X9 n) `