不等距离散点的积分程序
function = cubint(X,F,N,ia,ib)% 参考资料：Philio J. Davis 数值积分法 高等教育出版社 1986 p296-298
% X,F:积分点，函数值
% N: 积分点数，N>=4
% ia ib: ia,ib须为（1，N）整数 X(ia),X(ib)
% result: 积分值
% error: 误差估计
% ind: 条件满足为1，否则为0
% 对离散点三次插值，给出了基于四阶差商逼近四阶导数的误差估计
%
% 例：
% x=;
% y=x.^0.25;
% = cubint(x,y,16,1,16)；
% 运行结果：
% result = 0.7905 精确解：0.8
% error = -8.7208e-005
% ind = 1
ind=0;
if N<4 || ia<1 || ib>N
disp('wrong in N ia ib');
else
ind=1;
result=0;
error=0;
if ia==ib
disp('ia=ib,please check ia or ib');
else
if ia>ib
ind=-1;
it=ib;
ib=ia;
ia=it;
end
s=0;
c=0;
r4=0;
J=N-2;
if ia<(N-1) || N==4
J=max(3,ia);
end
K=4;
if ib>2 || N==4
K=min(N,ib+2)-1;
end
for i=J:K
if i<=J
H2=X(J-1)-X(J-2);
d3=(F(J-1)-F(J-2))/H2;
H3=X(J)-X(J-1);
d1=(F(J)-F(J-1))/H3;
H1=H2+H3;
d2=(d1-d3)/H1;
H4=X(J+1)-X(J);
r1=(F(J+1)-F(J))/H4;
r2=(r1-d1)/(H4+H3);
H1=H1+H4;
r3=(r2-d2)/H1;
if ia<=1
result=H2*(F(1)+H2*(0.5*d3-H2*(d2/6-(H2+H3+H3)*r3/12)));
s=-H2^3*(H2*(3*H2+5*H4)+10*H3*H1)/60;
end
else
H4=X(i+1)-X(i);
r1=(F(i+1)-F(i))/H4;
r4=H4+H3;
r2=(r1-d1)/r4;
r4=r4+H2;
r3=(r2-d2)/r4;
r4=(r3-d3)/(r4+H1);
end
if i>ib || i<=ia
error=error+r4*s;
else
result=result+H3*((F(i)+F(i-1))*0.5-H3*H3*(d2+r2+(H2-H4)*r3)/12);
c=H3^3*(2*H3^2+5*(H3*(H4+H2)+2*H2*H4))/120;
error=error+(c+s)*r4;
if i==j
s=s+c+c;
else
s=c;
end
end
if i>=K
if ib>=N
result=result+H4*(F(N)-H4*(0.5*r1+H4*(r2/6+(H3+H3+H4)*r3/12)));
error=error-H4^3*r4*(H4*(3*H4+5*H2)+10*H3*(H2+H3+H4))/60;
else
if ib>=N-1
error=error+s*r4;
end
end
else
H1=H2;
H2=H3;
H3=H4;
d1=r1;
d2=r2;
d3=r3;
end
end
if ind~=1
it=ib;
ib=ia;
ia=it;
result=-result;
error=-error;
ind=1;
end
end
end
end最近在做边界元,要积分
也看到有人问过离散点的积分
刚好,刚调完,结果正确
是根据fortran改过来的
抛砖引玉了
[ 本帖最后由 sigma665 于 2008-1-16 20:19 编辑 ]
页:
[1]