Labfans是一个针对大学生、工程师和科研工作者的技术社区。 论坛首页 | 联系我们(Contact Us)
MATLAB爱好者论坛-LabFans.com
返回   MATLAB爱好者论坛-LabFans.com > 工程数学软件 > MATLAB论坛
MATLAB论坛 一切MATLAB相关问题在此讨论。
回复
 
主题工具 显示模式
旧 2011-11-16, 10:36   #1
fanxiangn
初级会员
 
注册日期: 2011-10-14
年龄: 36
帖子: 2
声望力: 0
fanxiangn 正向着好的方向发展
默认 有个三次样条插值的程序不知道为何,误差有点大

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   #2
fanxiangn
初级会员
 
注册日期: 2011-10-14
年龄: 36
帖子: 2
声望力: 0
fanxiangn 正向着好的方向发展
默认 回复: 有个三次样条插值的程序不知道为何,误差有点大

中间追赶法的程序
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
fanxiangn 当前离线   回复时引用此帖
回复

主题工具
显示模式

发帖规则
不可以发表新主题
不可以发表回复
不可以上传附件
不可以编辑自己的帖子

启用 BB 代码
论坛启用 表情符号
论坛启用 [IMG] 代码
论坛禁用 HTML 代码



所有时间均为北京时间。现在的时间是 17:42


Powered by vBulletin
版权所有 ©2000 - 2025,Jelsoft Enterprises Ltd.