数值积分的计算

时间 2020-12-27

对于积分 $I=\int_{a}^{b}f(x)dx$ ,只要找到被积函数的原函数，便有下列牛顿--莱布尼兹公式：

$\int_{a}^{b}f(x)dx=F(b)-F(a)$

但实际上采用这种方法往往有困难，一是有些原函数不能用初等函数表达，有些计算过于麻烦。因为我们可以用数值的方法计算。

积分中值定理告诉我们，在区间 [a,b]内存在一点 $\delta$ ，使得

$\int_{a}^{b}f(x)dx=(b-a)f(\delta )$ .

就是说，底为b-a而高为 $f(\delta )$ 的矩形面积恰等于所求曲边梯形的面积 I ，问题在于点 $\delta$ 一般是不知道的。但我们可以将 $f(\delta )$ 称为区间上的平均高度，这样我们只要对平均高度提供一种算法，便可以获得一种数值求积方法。

一.小区间上的数值积分

1.梯形公式

我们找一个函数，即最简单的一次函数来拟合,用两端点 f(a),f(b) 的算术平均值作为平均高度的近似值，那么 $f(\delta )=(\frac{f(a)+f(b)}{2})$ ，那么求积公式为

$\int_{a}^{b}f(x)dx=\frac{b-a}{2}(f(a)+f(b))$

,这就是梯形公式。

代码如下：

%f(x)=2*x  
a=1;b=1.1;

disp("梯形公式解");
t=((b-a)/2)*(hanshu(a)+hanshu(b))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)
function[f]=hanshu(x)
  f=x*x;
end

结果：

梯形公式的代数精度为1，不再赘述。

2.辛普森公式

如果我们不用最简单的函数，而改用二次函数，那么便是辛普森公式，直接写公式

$\int_{a}^{b}f(x)dx=\frac{b-a}{6}(f(a)+f(b)+4f(\frac{a+b}{2}))$

disp("sinpson解")
t=((b-a)/6)*(hanshu(a)+hanshu(b)+4*hanshu((a+b)/2))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
  f=x*x;
end

结果：

同样的函数，辛普森公式的解就比梯形公式的解更精确了，当然计算复杂度也高了。

二：一般的数值积分

1. 复合的梯形公式

一般的区间主要是切分成小区间，再利用梯形公式，

对区间 [a,b]做切分，分成m块，Xi=(b-a)/m;

$\int_{a}^{b}f(x)dx=\sum_{1}^{m}\int_{x_{i-1}}^{x_{i}}f(x)dx$

分成m块：

$\sum_{1}^{m}\int_{x_{i-1}}^{x_{i}}f(x)dx=\sum_{1}^{m}(\frac{b-a}{2m})(f(x_{i-1})+f(x_{i}))$

利用梯形公式：（中间的边需要算2遍，所以乘2）

$\sum_{1}^{m}(\frac{b-a}{2m})(f(x_{i-1})+f(x_{i}))=\frac{b-a}{2m}(f(a)+f(b)+2*\sum_{1}^{m-1}f{(x_{i}})$

代码如下：

disp("复合的梯形公式,分成2块")
m=2;step=(b-a)/m;
t3=((b-a)/(2*m))*(hanshu(a)+hanshu(b)+xun(a,step,b))
disp("复合的梯形公式,分成5块")
m=5;step=(b-a)/m;
t3=((b-a)/(2*m))*(hanshu(a)+hanshu(b)+xun(a,step,b))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
  f=x*x;
end
function[t]=xun(a,step,b)
t=0;
  for i=a+step:step:b-step
      t=t+hanshu(i);
  end
  t=2*t;
end

结果：

2. 复合的辛普森公式

直接写公式：

$\int_{a}^{b}f(x)dx=\frac{b-a}{6m}(f(a)+f(b)+2*\sum_{1}^{m-1}f(x_{i})+4*\sum_{1}^{m}f(x_{_{i-\frac{1}{2}}}))$

代码如下：

disp("复合的sinpson公式,分成1块")
m=1;step=(b-a)/m;
t3=((b-a)/(6*m))*(hanshu(a)+hanshu(b)+xun(a,step,b)+xu(a,step,b))
disp("复合的sinpson公式,分成2块")
m=2;step=(b-a)/m;
t3=((b-a)/(6*m))*(hanshu(a)+hanshu(b)+xun(a,step,b)+xu(a,step,b))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
  f=x*x;
end
function[t]=xun(a,step,b)
t=0;
  for i=a+step:step:b-step
      t=t+hanshu(i);
  end
  t=2*t;
end
function[xx]=xu(a,step,b)
 xx=0;
   for i=a+step/2:step:b
       xx=xx+hanshu(i);
   end
   xx=xx*4;
end

结果如下：