计算方法大作业——龙贝格积分法的matlab程序

附录4 龙贝格积分法的matlab程序

clc;clear; s = input('请输入被积函数表达式：\\nf(x) = ','s'); a = input('\\n请输入积分区间下限：\\n'); b = input('\\n请输入积分区间上限：\\n'); e = input('\\n请输入积分误差：\\n'); N = -log10(e); f = inline(s); i = 1; error = 1; syms x f_a = limit(s,x,a); f_b = limit(s,x,b); T(1) = (b-a)/2*(double(f_a)+double(f_b)); S(1) = 0; C(1:2) = 0; R(1:3) = 0; while error>e || i<5 i = i+1; sigma = 0; for j=1:2^(i-2) sigma = sigma+f(a+(2*j-1)*(b-a)/2^(i-1)); end T(i) = 1/2*T(i-1)+(b-a)/2^(i-1)*sigma; S(i) = T(i)+1/3*(T(i)-T(i-1)); if i>=3 C(i) = S(i)+1/15*(S(i)-S(i-1)); if i>=4 R(i) = C(i)+1/63*(C(i)-C(i-1)); error = abs(R(i)-R(i-1)); end end end for k=1:i T_int = fix(T(k)); T_dec = round((T(k)-T_int)*10^N); S_int = fix(S(k)); S_dec = round((S(k)-S_int)*10^N); C_int = fix(C(k)); C_dec = round((C(k)-C_int)*10^N); R_int = fix(R(k)); R_dec = round((R(k)-R_int)*10^N); if k==1 fprintf('\\n%0.f.%0.f\\n',T_int,T_dec); elseif k==2 fprintf('%0.f.%0.f %0.f.%0.f\\n',T_int,T_dec,S_int,S_dec); elseif k==3 fprintf('%0.f.%0.f %0.f.%0.f %0.f.%0.f\\n',T_int,T_dec,S_int,S_dec,C_int,C_dec); else fprintf('%0.f.%0.f %0.f.%0.f %0.f.%0.f %0.f.%0.f\\n',T_int,T_dec,S_int,S_dec,C_int,C_dec,R_int,R_dec); end end fprintf('\\n所求积分值为 %0.f.%0.f,误差为%.e。\\n',R_int,R_dec,e);27