登录论坛

查看完整版本 : [MATLAB数学相关] 有个三次样条插值的程序不知道为何,误差有点大


fanxiangn
2011-11-16, 10:36
x=[8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2];
y=[0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177];
Y1=0.01087;
Y2=100;
n=size(x,2);
for i=1:n-2 %形成系数矩阵
u(i)=(x(i+1)-x(i))/(x(i+2)-x(i));
end
for i=1:n-2
la(i)=(x(i+2)-x(i+1))/(x(i+2)-x(i));
end
A=2*eye(n-2,n-2)+diag(u(1:n-3),1)+diag(la(2:n-2),-1);
for i=1:n-1 %形成非齐次项
f(i)=(y(i+1)-y(i))/(x(i+1)-x(i));
end
for i=1:n-2
g(i)=3*(la(i)*f(i)+u(i)*f(i+1));
end
g(1)=g(1)-la(1)*Y1;
g(n-2)=g(n-2)-u(n-2)*Y2;
disp(g);
[m,L,U]=zhuigan_Y(A,g);
m=[Y1;m';Y2];
syms x1 aerfa beita;
syms s s1;
for i=2:n
aerfa(i)=((x1-x(i-1))^2/(x(i)-x(i-1))^2)*(1+2*(x1-x(i))/(x(i-1)-x(i)));
end
aerfa(1)=(x1-x(2))^2/(x(1)-x(2))^2*(1+2*(x1-x(1))/(x(2)-x(1)));
for i=2:n
beita(i)=((x1-x(i-1))/(x(i)-x(i-1)))^2*(x1-x(i));
end
beita(1)=((x1-x(2))/(x(1)-x(2)))^2*(x1-x(1));
for i=1:n
s(i)=y(i)*aerfa(i)+m(i)*beita(i);
end
for i=1:n-1 %样条插值在各个区间上的表达式
s1(i)=s(i)+s(i+1);
end
fplot(inline(s1(2)),[x(2),x(3)]);
着个区间上的图形有问题,第一类边界条件在左边界一阶倒数为Y1,右边界一阶导数为Y2
希望老鸟们帮助帮助小弟

fanxiangn
2011-11-16, 10:38
中间追赶法的程序
unction [x,L,U] = zhuigan_Y(A,g)
% A:系数矩阵
% g:等号右端的向量
aerfa(1)=A(1,1); %产生追赶法的L,U矩阵
gama=diag(A,-1);
c=diag(A,1);
n=size(A,1);
beita=zeros(1,n-1);
for i=1:n-1
beita(i)=c(i)/aerfa(i);
aerfa(i+1)=A(i+1,i+1)-gama(i)*beita(i);
end
L=diag(aerfa)+diag(gama,-1);
U=eye(n,n)+diag(beita,1);
y=zeros(1,n);
y(1)=g(1)/A(1,1);
for i=2:n
y(i)=(g(i)-gama(i-1)*y(i-1))/(A(i,i)-gama(i-1)*beita(i-1));
end
x(n)=y(n);
for j=n-1:-1:1
x(j)=y(j)-beita(j)*x(j+1);
end
end