数学建模社区-数学中国
标题:
程序有点问题？？？？？？？帮帮忙
[打印本页]
作者:
潇湘夜雨123456
时间:
2012-12-15 11:22
标题:
程序有点问题？？？？？？？帮帮忙
function [f,f0] = ThrSample1(x,y,y_1,y_n,x0)
% v: }4 t7 ^( a
%求已知数据点的第一类三次样条差值多项式及其插值点处的值
+ p" ~1 G+ E; |6 G2 @
%已知数据点的x坐标向量：x
/ @- R6 a, X( p6 q. e$ y
%已知数据点的y坐标向量：y
5 j/ z ~1 G5 i' `' n/ l
%左端点的一阶导数：y_1
d( X$ r* E8 `# P
%右端点的一阶导数：y_n
1 W3 h3 y. L- Q+ O. u
%插值点的x坐标：x0
6 k+ A: S! B8 m! X% [
%求得的三次样条差值多项式：f
" d0 z: P# F. S
%求得的x0处的插值：f0
7 j# u# s& R4 t5 f% `; f
syms t
3 v7 @9 A* s/ }6 ?# a+ L
f = 0.0000;
0 e* R* J. ^+ y2 v# Y
f0 = 0.0000;
" D) r$ N6 M# H, ^* ~% O$ w
if(length(x)==length(y))
3 Y9 H* T! \; q/ N. }+ W
n = length(x);
/ q4 h; G$ s4 |9 }& j
else
2 e$ e; Z! i+ |2 c7 m6 L7 B8 e
disp('x和y的维数不相等！');
- h0 N; ], d% e2 _
return;
3 r2 v4 `; l1 P* X, k( s; K& [
end %维数检查
5 r1 ?+ t0 E$ E7 m, L
for i=1:n
+ c/ {% ?& z3 R$ z6 @
if(x(i)<=x0)&&(x(i+1)>=x0)
V0 y3 w* @0 z4 j9 O
ihdex = i;
7 l' r+ v+ q' }+ z. |; a# X: T2 N
break;
2 R. o' D k) p2 G$ F
end
% V2 a/ r, l+ F( i
end %找到x0所在区间
0 O. m2 E2 t4 |! q+ c+ I% u
A = diag(2*ones(1,n)); %求解m的系数矩阵
( U) f& x% D: f9 j% @
u = zeros(n-2,1);
2 k) L# X6 J) m& W
lamda = zeros(n-1,1);
) ` E& ~8 d4 r- j4 b8 E7 c5 n- X
c = zeros(n,1);
) {- ?& u0 v$ J/ m; y4 |
for i = 2:n-1
$ a$ m% {+ i0 t; v
u(i-1) = (x(i)-x(x-1))/(x(i+1)-x(i-1));
5 K! N' j( Z$ Y" T% ]% H5 f
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
' ?, B& I5 T/ S2 i
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));
2 d: m: F3 t' \, u
A(i,i+1) = u(i-1);
( @( d/ }# ?7 U4 m" L- |1 j* q
A(i,i-1) = lamda(i); %构造系数矩阵及向量c
, q3 _) \; V+ H8 h! {, i% p
end
" e2 X! ]) @8 ~* ]$ S8 X0 T
c(1) = 2*y_1;
. U+ s6 V0 e$ y( z8 x" o
c(n) = 2*y_n;
# Q/ W0 [6 T$ k) U4 X4 m5 z* x7 w
m = followup(A,c); %用追赶法求解方程组
* o& b2 C8 y* \# _# V& l# M
h = x(index+1) - x(index); %x0所在区间长度
; s- S. |# f, J
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...
3 t# B3 f8 o; @/ C
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...
, ^! j: |# q) m- ?1 ]' q
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...
) i( B) `3 K/ J o. y
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h; %x0所在区间的插值函数
" R+ K4 C- a4 A% K4 `- @
f0 = abs(f,'t',x0) %x0处的插值
1 d4 g# L% @, U, U
# w9 y1 j; \$ U
, @& b. L6 a. m; g" `
' ?* }8 J4 T: W4 o
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5