数值积分方法：复合梯形公式、复合辛普森公式和龙贝格算法

使用MATLAB代码实现数值积分方法

本文将使用MATLAB代码演示三种常见的数值积分方法：复合梯形公式、复合辛普森公式和龙贝格算法。

1. 定义函数

function y=fun41(x)
y=sqrt(x).*log(x);
end

该代码定义了一个名为fun41的函数，输入为x，输出为sqrt(x)乘以log(x)的值。

2. 复合梯形公式和复合辛普森公式

clear;
clc;
h=0.001; %h为步长,可分别令h=1,0.1,0.01,0.001
n=1/h;
t=0;
for i=1:n-1
  t=t+fun41(i*h);
end
T=h/2*(0+2*t+fun41(1));
T=vpa(T,10) % 以上为复合梯形公式

% 以下为复合辛普森公式
s1=0;
s2=0;
for i=0:n-1
  s1=s1+fun41(h/2+i*h);
end
for i=1:n-1
  s2=s2+fun41(i*h);
end
S=h/6*(0+4*s1+2*s2+fun41(1));
S=vpa(S,10)

该代码段首先使用复合梯形公式计算数值积分，然后使用复合辛普森公式计算数值积分。代码解释如下：

clear; clc;：清空命令窗口和工作空间。
h=0.001; n=1/h; t=0;：设定步长h为0.001，n为1/h，t为0。
for i=1:n-1 t=t+fun41(i*h); end：循环n-1次，每次将fun41(i*h)加到t中。
T=h/2*(0+2*t+fun41(1)); T=vpa(T,10)：使用复合梯形公式计算数值积分T，公式为h/2*(0+2*t+fun41(1))。使用vpa函数将T保留10位有效数字。
s1=0; s2=0;：初始化变量s1和s2为0。
for i=0:n-1 s1=s1+fun41(h/2+i*h); end：循环n-1次，计算s1的值。
for i=1:n-1 s2=s2+fun41(i*h); end：循环n-1次，计算s2的值。
S=h/6*(0+4*s1+2*s2+fun41(1)); S=vpa(S,10)：使用复合辛普森公式计算数值积分S，公式为h/6*(0+4*s1+2*s2+fun41(1))。使用vpa函数将S保留10位有效数字。

3. 龙贝格算法

clear;
clc;
m=16;
h=1;
T(1)=(0+fun41(1))*h/2
for i=2:m
  h=h/2;
  n=1/h;
  t=0;
  for j=1:2:n-1
    t=t+fun41(j*h);
  end
  T(i)=T(i-1)/2+h*t; %梯形公式
end
for i=1:m-1
  for j=m:i+1
    T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1);
    %通过不断的迭代求得T(j)，即T表的对角线上的元素。
  end
end
vpa(T(m),10)

该代码段使用龙贝格算法计算数值积分。代码解释如下：

clear; clc;：清空命令窗口和工作空间。
m=16; h=1; T(1)=(0+fun41(1))*h/2：设定m为16，h为1，将T表的第一个元素设置为(0+fun41(1))*h/2。
for i=2:m h=h/2; n=1/h; t=0; for j=1:2:n-1 t=t+fun41(j*h); end T(i)=T(i-1)/2+h*t; end：循环m-1次，每次将h减半，n设为1/h，t设为0，使用梯形公式计算当前T表中的元素，公式为T(i)=T(i-1)/2+h*t。
for i=1:m-1 for j=m:i+1 T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1); end end：循环m-1次，使用龙贝格算法迭代求解T表中的对角线元素，公式为T(j)=4^i/(4^i-1)*T(j)-1/(4^i-1)*T(j-1)。
vpa(T(m),10)：使用vpa函数将T表的最后一个元素保留10位有效数字。

注意： 本文仅提供示例代码，具体的代码实现可能需要根据实际情况进行修改。

希望本文能够帮助您理解数值积分方法的实现原理。如果您有任何疑问，请随时提出。