用三次樣條插值求斜率
三次樣條插值的matlabt程序function x=followup(a,b,c,d)n=length(d);5 @2 m s5 {8 ^: C$ L
a(1)=0;
%“追”的過程2 o) O( M: o& Z' ]
L(1)=b(1);$ n3 }2 S: F* c
y(1)=d(1)/L(1);
u(1)=c(1)/L(1);0 a; _" y. [0 x3 d6 l
for i=2
n-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);
end! n+ w" {8 ^ V
L(n)=b(n)-a(n)*u(n-1);
y(n)=(d(n)-y(n-1)*a(n))/L(n);
%“趕”的過程
x(n)=y(n);- [- Q' O. v/ i9 F1 O
for i=(n-1):-1:1
x(i)=y(i)-u(i)*x(i+1);, S. A8 }; N: m* t% X
end
5 ^. O( F: l* w2 a
function[s,y0]=spline3 (x,y,x0)( _+ v9 g# C; s" w+ I6 f; G+ \, {
%x,y為數表x0為插值點s表示插值函數y0為x0對應的插值函數值2 ~2 H/ O" q9 }: P6 d4 O$ e
syms t, _7 _: v5 a! b4 c& Z
n=length(x);/ I6 D# P! o' C9 c5 G5 `! a
%得出n
for i=1:n-1;
h(i)=x(i+1)-x(i);+ u: a, x* ]3 h, k
end
for i=2:n-1;
lamda(i)=h(i)/(h(i-1)+h(i));) V1 s* c6 U$ U0 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)));+ ~ ~9 y+ E1 g# j: F3 m' t8 `
end
g(1)=3*((y(2)-y(1))/h(1));
g(n)=3*((y(n)-y(n-1))/h(n-1));
%前邊求出lamda miu和g從而可以確定系數矩陣1 x' X0 @. S4 Z1 x
miu(1)=1;
miu(4)=0;* j" M* p; j1 P
lamda(n)=1;
lamda(1)=0;
%根據第二邊界條件補充兩個lamda和miu的值! f7 w* S+ _0 l/ H: m6 }
for i=1:n
beta(i)=2;1 t& n( o4 S! I
end6 h- b/ W" W+ k$ c
m=followup(lamda,beta,miu,g)/ {/ H" r) K$ e# Z
%解出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)...# ]4 S" W0 @- i+ V! H
+(t-x(i))^2*(h(i)+2*(x(i+1)-t))*y(i+1)/(h(i)^3)...2 y0 I+ S5 L! x; R
+(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 d2 E0 ?2 y! L# G1 g7 H
end/ S8 [! r* p- V: d2 Q
%得到插值的結果各段的t的表達式
%接下來要將插值點x0代入首先確定x0所在的插值區間' z; O/ h/ s) E" j' q
for i=1:n-1
if (x(i)>x0)9 r6 z7 K* E6 z7 N2 s% ~
in=i;- s I3 @1 x) ?5 z3 b6 E
end
end
s=st(in);
s=expand(s);
s=collect(s,'t');
y0=subs(s,'t',x0)
%s是插值多項式y0是插值點的函數值0 W7 q- N7 c- B! B. Y( e
在matlab中輸入
x=[1 2 4 5]；
y=[1 3 4 2]；
spline3(x,y,2)
會得到各點的斜率
# x% g7 F4 Z$ X: a# i
" _" r4 F" }% ~3 B' `