第一部分:问题分析
(1)实验题目:龙贝格积分算法
具体实验要求:用matlab编写龙贝格积分的代码,要求代码实现用户输入了被积函数、积分区间、精度之后,龙贝格积分表(T-数表)。
(2)实验目的:让同学们进一步掌握龙贝格积分的原理以及运算过程,并且通过matlab编程培养实际的上机操作能力和代码能力。
第二部分:数学原理
龙贝格求积公式也称为逐次分半加速法。它是在梯形公式、辛普森公式和柯特斯公式之间的关系的基础上,构造出一种加速计算积分的方法。 作为一种外推算法,它在不增加计算量的前提下提高了误差的精度。
在等距基点的情况下,用计算机计算积分值通常都采用把区间逐次分半的方法进行。这样,前一次分割得到的函数值在分半以后仍可被利用,且易于编程。
根据已知的龙贝格外推公式:
第三部分:程序设计流程
- 功能函数(被调用者)
注:1.代码特别之处,将积分表作为一个矩阵存储,在打印时直接打印
2.要显示数据的精度,需要读者在命令行中输入format long指令
注意:该算法的T数表和实验原理中的T数表的顺序一样,但计算过程是一样的,只是修改了排列顺序。(一定不要把龙贝格积分算法中,代表各步骤的积分式T(m,k)和矩阵元素A(i,j)混淆)
- 主函数(面向用户:调用者)
第四部分:代码实现
Romberg积分函数实现:
function [] = Romberg_Iteration(f,a,b,e)
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a;
T=double(h/2*(f(a)+f(b)));%梯形公式求出T(1,1)
err=b-a;
while err>=e
k=k+1;
h=h/2;
tmp=0;
for i=1:n
tmp=tmp+f(a+(2*i-1)*h);
end
T(k+1,1)=double(T(k)/2+h*tmp);%求出行首元
for j=1:k
T(k+1,j+1)=double(T(k+1,j))+double((T(k+1,j)-T(k,j))/(4^j-1));%迭代算法
end
n=n*2;
err=abs(T(k+1,k+1)-T(k,k));%误差为该次迭代的首元和上一次迭代首元的差
end
disp(T);
用户调用部分:
用户在该函数函数中输入积分区间的左右端点a、b,精度e,以及被积函数f(x)
a = 0;
b = 1;
e = 1e-6;
syms x;
f(x) = exp(x)*x^2;
Romberg_Iteration(f,a,b,e);