2 详解matlab实现龙贝格算法求复杂被积函数的积分

目录

2 习题二

2.1 题目

2.2 问题背景

2.3 算法

2.4 matlab编程实现

2.5 结果展示

2.6 讨论

---

2 习题二

2.1 题目

编写Romberg算法通用程序计算下列积分，允许误差10-6

2.2 问题背景

实际计算时遇到的被积函数往往很复杂，找不到相应的原函数；即使是一些简单的函数，比如e-x2 都找不到用初等函数表示的原函数。另外，在一些计算问题中，f(x) 的值是通过观测或数值计算得到一组数据表，此时显然不能用 Newton-Leibniz 公式计算积分。所有这些因素，促进人们去研究定积分的近似计算，定积分的近似计算称为数值积分。构造数值积分公式最通常的方法是用积分区间上的 n 次插值多项式代替被积函数，由此导出的求积公式称为插值型求积公式。例如梯形公式 (TrapezoidalApproximations) 与抛物线公式 (Approximations Using Parabolas) 就是最基本的近似公式。但它们的精度较差。

龙贝格算法是在区间逐次分半过程中，对梯形公式的近似值进行加权平均获得准确程度较高的积分近似值的一种方法，它具有公式简练、计算结果准确、使用方便、稳定性好等优点，因此在等距情形宜采用龙贝格求积公式 (Rhomberg Integration)

2.3 算法

2.4 matlab编程实现

（将函数文件放在与主程序同一个文件夹下，然后运行主程序）

主程序：

clear;

format long

fun1=@(x)x*sin(x);   %第一问的被积函数

a1=0;b1=2*pi;  %第一问的积分区间

fun2=@(x)2/sqrt(pi)*exp(-x^2);   %第二问的被积函数

a2=0;b2=1;  %第二问的积分区间

eps=10^(-6);  %允许误差

[I1,R1]=romberg(fun1,a1,b1,eps)  %I1为第一问的积分近似值，R1为存储所计算第一问Rn值的矩阵

[I2,R2]=romberg(fun2,a2,b2,eps)  %I2为第二问的积分近似值，R2为存储所计算第二问Rn值的矩阵

%[I3,R3]=romberg(@(x)exp(-x^2),0,1,0.5*10^(-5)) ;%课本例题验证

函数：

function [I,R]=romberg(fun,a,b,eps)

%龙贝格积分法求积分近似值

%fun表示函数表达式，a是积分下限，b是积分上限，esp是允许误差

format long

syms f(x)

f(x)=fun;

T=zeros(10,1);S=zeros(10,1);C=zeros(9,1);R=zeros(8,1);

T(1)=((b-a)/2)*(f(b)+f(a));

for i=1:5

  %计算前五行T_1到T_32

  e=0;

  for j=1:2^(i-1)

      e=e+f(a+(2*j-1)*(b-a)/(2^i));

  end

  T(i+1)=(T(i)+(2^(1-i)*(b-a))*e)/2;

end

for k=1:4

    S(k)=(4*T(k+1)-T(k))/3;

end

for m=1:3

   C(m)=(16*S(m+1)-S(m))/15;

end

for n=1:2

    R(n)=(64*C(n+1)-C(n))/63;

end

while abs(R(2)-R(1))<eps  %第一次判断迭代终止条件|R_2-R_1|<eps

   I=R(2);

   break;

end

%若R2不满足条件，再继续往下算

for i=6:10

  e=0;

  for j=1:2^(i-1)

      e=e+f(a+(2*j-1)*(b-a)/(2^i));

  end

  T(i+1)=(T(i)+(2^(1-i)*(b-a))*e)/2;

end

for k=5:10

    S(k)=(4*T(k+1)-T(k))/3;

end

for m=1:9

   C(m)=(16*S(m+1)-S(m))/15;

end

for n=3:8

    R(n)=(64*C(n+1)-C(n))/63;

    if (n>1)&&(abs(R(n)-R(n-1))<eps)  %迭代终止条件|R_2n-R_n|<eps

        I=R(n);

        break;

    end

end

end

2.5 结果展示

>> work2

I1 =

  -6.283185307367059

R1 =

  -6.266954014124826

  -6.283202742954948

  -6.283185358349977

  -6.283185307367059

                   0

                   0

                   0

                   0

I2 =

   0.842700792949092

R2 =

   0.842700663941961

   0.842700792763488

   0.842700792949092

                   0

                   0

                   0

                   0

                   0

2.6 讨论

（1）

这里将Rn的计算分为两部分，以降低时间复杂度。因为通常情况下R2的精度已经很高，能够满足要求的情况下就不用继续往下算。

（2）

例题验证调用代码，在主程序末尾添加以下代码

[I3,R3]=romberg(@(x)exp(-x^2),0,1,0.5*10^(-5)) ;

结果显示误差数量级为1*10^(-7)  ，误差极小，而允许误差为0.5*10^(-5)，误差可忽略不计。