数学建模社区-数学中国
标题:
程序有点问题？？？？？？？帮帮忙
[打印本页]
作者:
潇湘夜雨123456
时间:
2012-12-15 11:22
标题:
程序有点问题？？？？？？？帮帮忙
function [f,f0] = ThrSample1(x,y,y_1,y_n,x0)
+ a @ X6 c' G6 t3 J% i
%求已知数据点的第一类三次样条差值多项式及其插值点处的值
" Y- b# y9 i, `' L* k
%已知数据点的x坐标向量：x
/ p- j. _; o+ `) L4 q
%已知数据点的y坐标向量：y
+ h( F) w9 l7 b( _; O0 s1 a: I
%左端点的一阶导数：y_1
2 j( k. L+ N( v7 M/ z
%右端点的一阶导数：y_n
( s7 x( W+ o: Q7 A3 H
%插值点的x坐标：x0
; d( Z" r( ]3 J- e( @9 R. C% [
%求得的三次样条差值多项式：f
0 j; v; m' y. S3 T: `; d
%求得的x0处的插值：f0
2 v0 d! t* |! S2 V, Q" K3 J8 |
syms t
7 G, D* ^, O! p0 Y5 u$ l" D' ]' _
f = 0.0000;
. y1 V6 S" \9 b9 @% {
f0 = 0.0000;
) T& K7 G& j5 m3 l0 {; P
if(length(x)==length(y))
9 N% _! o4 ~' a& d& K4 ]
n = length(x);
" y. j' _1 p" ~. k! \
else
4 v3 K; P% k# B+ h/ | {
disp('x和y的维数不相等！');
) x. Q8 ^& h! _& N( }6 i
return;
6 v6 Z6 m( g5 |# R
end %维数检查
i7 @, W* n) I
for i=1:n
+ o3 k' B7 k1 S+ W: N( e5 U; G
if(x(i)<=x0)&&(x(i+1)>=x0)
/ n. S# C3 j6 i7 N) i! {
ihdex = i;
' ^ V D! P2 _) Z9 R
break;
2 T. z2 L& s3 h9 R3 M) X( H8 X, M
end
, L# [% e, u; M
end %找到x0所在区间
' r2 j) n r) s' C' ~( E! I% r
A = diag(2*ones(1,n)); %求解m的系数矩阵
5 J$ l3 M; C- T$ y
u = zeros(n-2,1);
$ E0 t1 Y& r! J6 f; n
lamda = zeros(n-1,1);
$ E7 f8 S' H) ]( K, V) b/ U3 C
c = zeros(n,1);
! B- I# n- G) a& Y
for i = 2:n-1
* R9 U9 ^ e3 S
u(i-1) = (x(i)-x(x-1))/(x(i+1)-x(i-1));
. h1 _2 o4 x0 F+ S4 ?
lamda(i) = (x(i+1)-x(i))/(x(i+1)-x(i-1));
, @& t+ @5 f& D l# H0 ]* Y
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));
V0 D- [2 {/ U/ u; k
A(i,i+1) = u(i-1);
) d% Y$ I/ n- t
A(i,i-1) = lamda(i); %构造系数矩阵及向量c
9 |9 y3 ?4 ]- A6 K5 J- ^; O
end
+ A, o% `$ q/ b& j
c(1) = 2*y_1;
$ ?0 U) F9 q6 J: C1 k
c(n) = 2*y_n;
& s' l, m' {- g
m = followup(A,c); %用追赶法求解方程组
. C0 W- Q1 J+ H5 `/ s; E) W6 i
h = x(index+1) - x(index); %x0所在区间长度
/ H) u; [: t$ V, }9 b
f = y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+...
6 k5 T9 O* w$ A# ^/ t. o+ c1 b
y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+...
( g+ d( O7 x0 F' A7 G" \ e
m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-...
, \5 r' l% p3 u/ ?% h4 A* S
m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h; %x0所在区间的插值函数
0 T3 E! u. j- c! ~. o$ v6 A
f0 = abs(f,'t',x0) %x0处的插值
2 F6 X) i# h* g- }
5 {- E7 J0 {. U( Y- S
/ J3 o9 a% Z3 ^" m
- o# p1 d0 U5 E/ c @
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5