数学建模社区-数学中国
标题:
程序有点问题？？？？？？？帮帮忙
[打印本页]
作者:
潇湘夜雨123456
时间:
2012-12-15 11:22
标题:
程序有点问题？？？？？？？帮帮忙
function [f,f0] = ThrSample1(x,y,y_1,y_n,x0)
" q: J' i* T. t; b* q
%求已知数据点的第一类三次样条差值多项式及其插值点处的值
' V- R5 H0 @1 }8 e2 B& _
%已知数据点的x坐标向量：x
. w/ c" n; g: V; |
%已知数据点的y坐标向量：y
, x1 s9 z1 F9 d0 R
%左端点的一阶导数：y_1
4 o9 r, _4 k H$ k6 I$ h# L1 X
%右端点的一阶导数：y_n
; m7 W. D* G$ W! G% S7 \
%插值点的x坐标：x0
& }' L0 W5 `, s' J% ^8 Q9 e
%求得的三次样条差值多项式：f
/ e: E. R5 `4 z
%求得的x0处的插值：f0
' ~0 o" V) p: U( m/ p3 P
syms t
7 P; C; {+ T' B8 l% @8 M- ~
f = 0.0000;
- _: Q% [1 E. L4 E' d. K3 Q9 V: w, \
f0 = 0.0000;
, W3 H( K/ r: g7 Z. R7 B1 k
if(length(x)==length(y))
% U8 i5 {; S2 t/ r* T# S1 v- o
n = length(x);
! R/ t1 B; w$ v: K: \% ~: x4 x- W
else
2 p, Z* h/ E) H- N
disp('x和y的维数不相等！');
. q& f( b7 ~! n) q- B
return;
- ^9 b' g0 P5 q3 b7 L) R7 D
end %维数检查
* k% F' N; x% g' @: ]8 ]6 ]
for i=1:n
v$ y0 y7 d7 ?: Z
if(x(i)<=x0)&&(x(i+1)>=x0)
1 ~/ S3 X) B, Z
ihdex = i;
: o, ]4 Q0 ^$ t3 @% o
break;
$ r. J( ?0 L, B! B% V
end
9 ^. W0 b( c. C3 {+ J
end %找到x0所在区间
0 e: O# V7 W) P2 I, u" c
A = diag(2*ones(1,n)); %求解m的系数矩阵
# h' q0 G5 N5 \7 X; d3 J' d
u = zeros(n-2,1);
- E2 y4 w% N5 Z+ x" p
lamda = zeros(n-1,1);
; p- E- H+ }/ l4 R: C. Y4 N
c = zeros(n,1);
! T+ i7 o/ j7 n! [. E( j
for i = 2:n-1
2 q. A1 d: x0 k" d) }
u(i-1) = (x(i)-x(x-1))/(x(i+1)-x(i-1));
. x" f+ O. ]; C, _7 L
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
\2 r( d6 v) w7 m* t0 [0 f
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/ L8 w" R+ g [# w4 `" i- J5 m3 V$ h
A(i,i+1) = u(i-1);
: j+ v# p& X1 l9 Z2 U* F
A(i,i-1) = lamda(i); %构造系数矩阵及向量c
; z( D* c3 D: h; {, K4 A) u+ V' V
end
. r4 B# Q' v" f& h- s! d7 v
c(1) = 2*y_1;
7 Y; F6 C' f+ W
c(n) = 2*y_n;
. j3 v8 B, G, _! t2 ]) Y# n, D
m = followup(A,c); %用追赶法求解方程组
! e% w+ W) @% `$ `' o9 [
h = x(index+1) - x(index); %x0所在区间长度
# ~6 D: j* U+ h
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...
( e* n' M; `+ Y" h6 {3 ?
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...
$ O' f) h7 C v/ Y# @6 {8 f
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...
/ E2 }; E) z# N" p3 [
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h; %x0所在区间的插值函数
! E1 q, y m: v$ ~
f0 = abs(f,'t',x0) %x0处的插值
. U# f2 r5 N2 s
; p& l( b, K; |. g8 ?9 @5 D2 s
. |& X4 C$ t% c1 n- N/ g
, D9 U% M* Q# {; \# ~/ ~
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5