数学建模社区-数学中国
标题:
程序有点问题？？？？？？？帮帮忙
[打印本页]
作者:
潇湘夜雨123456
时间:
2012-12-15 11:22
标题:
程序有点问题？？？？？？？帮帮忙
function [f,f0] = ThrSample1(x,y,y_1,y_n,x0)
R: U& W# a- k( v2 P6 e
%求已知数据点的第一类三次样条差值多项式及其插值点处的值
1 _" t9 V- Q" A; Y
%已知数据点的x坐标向量：x
4 x" R+ G" e1 g
%已知数据点的y坐标向量：y
1 x3 H3 M( ]3 o# \8 E
%左端点的一阶导数：y_1
9 n, g2 v2 S3 C3 C) l; Z
%右端点的一阶导数：y_n
& ]% Y7 N" ~8 |1 y- b( H5 V
%插值点的x坐标：x0
1 H" z( M! x6 W8 M E$ e" s( p
%求得的三次样条差值多项式：f
* o8 X. V) M0 f7 \1 }7 [( o
%求得的x0处的插值：f0
0 f4 I) I# b' o" x* @
syms t
! l% h" o! g' A; J+ u, y3 u
f = 0.0000;
+ g9 }9 Z9 {. |) Y% M* Z' g% O
f0 = 0.0000;
. }2 H, [: T6 d7 |, j% k
if(length(x)==length(y))
) _; M# y, Q0 u+ |+ n( m
n = length(x);
9 B u2 U+ w8 Y3 M7 }
else
9 F; o- k$ e: A2 k# S
disp('x和y的维数不相等！');
& Y2 ?' |5 e+ v1 K+ I- V7 A# h
return;
( V) w9 m2 o) S6 W2 _& J7 e# F
end %维数检查
0 p, ^) Q1 K" `; W6 c
for i=1:n
! E- c* \( w0 m8 u6 I
if(x(i)<=x0)&&(x(i+1)>=x0)
$ E4 O6 w: G, e" | T; @5 B: M
ihdex = i;
% W7 S* E, H0 h: w6 y$ B
break;
% p! X2 E+ M! ^* i9 t1 R. n4 ^# L- U3 g9 g
end
$ z" H' D# k1 d
end %找到x0所在区间
- K& z: M. o, k2 f+ Y
A = diag(2*ones(1,n)); %求解m的系数矩阵
* j2 _5 z. S0 e% U
u = zeros(n-2,1);
+ M4 i) K# v; M3 W- H
lamda = zeros(n-1,1);
* W! N% ]2 m& G; G: d
c = zeros(n,1);
, \, v4 A Z, {3 y) G: u' K. ?" D
for i = 2:n-1
) Z0 p7 X$ f, _) }! u. k
u(i-1) = (x(i)-x(x-1))/(x(i+1)-x(i-1));
, A: h4 v# ~" L
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
+ a$ |, c# C/ a. w/ \
c(i) = 3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));
* \. V) I! l" m2 J3 \: e
A(i,i+1) = u(i-1);
6 I' }& c) z5 H. B% g
A(i,i-1) = lamda(i); %构造系数矩阵及向量c
; N& M7 u$ x8 g
end
9 ^6 O* n. x6 x1 F0 `
c(1) = 2*y_1;
( n8 m3 k* h$ R8 e" ?, w3 `9 d0 i
c(n) = 2*y_n;
& W/ t) E8 e+ W
m = followup(A,c); %用追赶法求解方程组
: Z2 D( A! M7 M
h = x(index+1) - x(index); %x0所在区间长度
0 z. V0 r, `: M$ H- `
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...
/ x; B7 E- F( N$ o+ f% Q. p
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...
9 @9 m' f; Z/ V
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...
) S# o/ u2 V. r, S! ^6 Z
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h; %x0所在区间的插值函数
9 `, e0 _) D& W3 S
f0 = abs(f,'t',x0) %x0处的插值
; G- C3 F. i" p. ^1 Q5 z
$ z3 u7 s; [( ^4 y5 k4 z
, A8 n0 J @. {+ \
& ?: a# N6 ^ h. ~* p
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5