数值计算day8-数值积分

上节课主要介绍了计算微分的几种数值方法，对一阶微分，最简单的莫过于两点前向差分、后向差分和中心差分这三种方法，其中中心差分的精度最高，这三种差分公式都可以通过推导泰勒展开式得到，而通过泰勒展开式还可以推导出三点前向差分和三点后向差分。对二阶微分，则可以推导出三点中心差分、三点前向差分、三点后向差分公式。这些数值方法均可以推广到数值偏微分的领域，且对于两个精度不高的数值算法，还可以使用Richardson外推加速算法得到一个精度更高的算法。本节课主要介绍计算积分的数值方法。

1. 矩形和中点法

区间([a,b])上的积分表示的是曲线(f(x))在([a,b])内的面积，若将([a,b])分解为(n)个小区间：(a=x_0<x_1<cdots<x_n=b)，(J_i=[x_{i-1},x_i])区间内的面积可以估算为一个黎曼和：(A_i=f(x^*_i)Delta x_i)，其中(x^*_i)表示(J_i)内的某个点，(Delta x_i = x_i-x_{i-1})，此时积分即为黎曼求和的极限：$$ I(f)=int^b_{a}f(x)dx=lim_{n ightarrow infty,Delta x ightarrow 0}sum^n_{i=1}f(x*_i)Delta x_i $$ 其中(Delta x = max {Delta x_i})

1.1 左点法与右点法

左点法是一种矩形法，每个小区间([x_{i-1},x_i])的面积使用矩形公式来计算，其中高取为区间左端点的函数值(f(x_{i-1}))，宽取为小区间的长度(Delta x_i = x_i - x_{i-1})，(Delta x_i = x_i-x_{i-1}), $$ I_{i}(f)=f(x_{i-1})Delta x_i $$

[I(f)=sum^n_{i=1}I_i(f)=sum^n_{i=1}f(x_{i-1})Delta x_i ]

右点法与左点法相对，高取为右端点的函数值(f(x_i))，

[I_{i}(f)=f(x_{i})Delta x_i ]

[I(f)=sum^n_{i=1}I_i(f)=sum^n_{i=1}f(x_{i})Delta x_i ]

1.2 中点法

中点法使用区间([x_{i-1},x_i])的中点函数值(f(frac{x_{i-1}+x_i}{2}))作为高，区间的面积为：$$I_i(f)=f(frac{x_{i-1}+x_i}{2})Delta x_i$$

[I(f)=sum^n_{i=1}I_i(f)=sum^n_{i=1}f(frac{x_{i-1}+x_i}{2})Delta x_i ]

下面以线性函数(y=x)为例，讨论如何控制左点法的结果在想要的误差范围内。区间([a,b])取为([0,1])，(y=x)在该区间上积分的真实值为(frac{1}{2})，利用左点法计算，将([0,1])等分为(n)个小区间，每个小区间的宽度为(h=frac{1}{n})，左点法计算的积分值为：$$h*(frac{1}{n}+frac{2}{n}+...+frac{n-1}{n})=frac{n-1}{2n}$$ 计算误差为(frac{1}{2n})，如果想要将误差的量级控制在(10^{-4})以内，则 至少需要将区间分为(10^4)个小区间。而在该例中，中点法计算是没有误差的。

2. 梯形法（Trapezoidal Method）

梯形法是对矩形法的一种改善，采用线性函数来计算积分，对点((a,f(a)))和((b,f(b)))进行牛顿插值，得到：$$f(x)=f(a)+(x-a)f[a,b]=f(a)+(x-a)frac{f(b)-f(a)}{b-a}$$ 将区间([a,b])的积分值计算为梯形面积，有：$$I(f)=f(a)(b-a)+frac{1}{2}f(b)-f(a)=frac{f(a)+f(b)}{2}(b-a)$$ 对区间进行划分(a=x_0<x_1<x_2<cdots<x_n=b)，则：$$I_i(f)=int_{x_{i-1}}^{x_i}f(x)dx=frac{f(x_{i-1})+f(x_i)}{2}(x_i-x_{i-1})$$ $$I(f)=sum^{n_{i=1}I_i(f)=sum}n_{i=1}frac{f(x_{i-1})+f(x_i)}{2}(x_i-x_{i-1})=frac{h}{2}[f(a)+f(b)]+hsum^{n-1}_{i=2}f(x_i)$$

MATLAB实现：

function I = trapezoidal(Fun,a,b,N)
% trapezoidal numerically integrate using the Composite Trapezoidal Method.
% Input Variables:
% Fun Name for the function to be integrated.
% (Fun is assumed to be written with element-by-element calculations.).
% a  Lower limit of integration.
% b  Upper limit of integration.
% N  Number of subintervals.
% Output Variable:
% I  Value of the integral.
 
h = (b-a)/N;
x=a:h:b;
F=Fun(x);
I=h*(F(1)+F(N+1))/2+h*sum(F(2:N));

3. Simpson 方法

3.1 Simpson 1/3 法

梯形法是采用线性函数插值，来估算区间面积，一种改进方法是使用非线性插值，包括二次插值、三次插值。Simpson 1/3 法采用二次插值法。

取区间([a,b])内的三点(x_1=a,x_3=b,x_2=frac{a+b}{2})，使用牛顿多项式插值法，过这三个点的曲线方程为：$$p(x)=alpha+eta (x-x_1)+gamma (x-x_1)(x-x_2)$$ 求解可得 $$alpha=f(x_1),eta = frac{f(x_2)-f(x_1)}{x_2-x_1}=frac{f(x_2)-f(x_1)}{h}$$,$$gamma = frac{f[x_3,x_2]-f[x_2,x_1]}{x_3-x_1}= frac{frac{f(x_3)-f(x_2)}{h}-frac{f(x_2)-f(x_1)}{h}}{2h}=frac{f(x_3)-2f(x_2)+f(x_1)}{2h^2}$$ 使用曲线(p(x))代替原曲线计算积分，$$I(f)=alpha(x_3-x_1)+eta(x_3-x_1)^2*frac{1}{2}+gamma int^{x_3}{x_1}(x-x_1)(x-x_2)dx$$ $$=2hf(x_1)+frac{f(x_2)-f(x_1)}{h}frac{4h^2}{2}+...=frac{h}{3}[f(x_1)+4f(x_2)+f(x_3)] $$
现在将区间([a,b])等分成(N)(偶数)个小区间：$$a=x_1<x_2<x_3<cdots < x_N<x{N+1}=b$$，则区间([x_{i-1},x_{i+1}])内的积分值为：$$ I_i(f)=frac{h}{3}[f(x_{i-1})+4f(x_i)+f(x_{i+1})]$$

对所有小区间积分值求和，有：$$ I(f)=frac{h}{3}[f(a)+4f(x_2)+f(x_3)+f(x_3)+4f(x_4)+f(x_5)+...+f(x_{N-1})+4f(x_N)+f(b)] $$

[=frac{h}{3}[f(a)+4sum_{i=2,4,6}^Nf(x_i)+2sum_{j=3,5,7}^{N-1}+f(b)] ]

3.2 Simpson 3/8 法

在Simpson 3/8 法中，对每四个点进行三次插值：$$p(x)=c_3x^3+c_2x2+c_1x+c_0$$ 小区间的积分值为：$$frac{3}{8}h[f(a)+3f(x_2)+3f(x_3)+f(b)]$$ 因此将区间([a,b])等分为(N)(3的倍数)个小区间，对每个小区间进行三次插值计算积分后求和，可得：$$ I(f)=frac{3h}{8}[f(a)+3f(x_2)+3f(x_3)+f(x_4)+f(x_4)+3f(x_5)+3f(x_6)+f(x_7)+ $$

[...+f(x_{N-2})+3f(x_{N-1})+3f(x_N)+f(b)] ]

[=frac{3h}{8}[f(a)+3sum^{N-1}_{i=2,5,8}[f(x_i)+f(x_{i+1})]+2sum^{N-2}_{j=4,7,10}f(x_j)+f(b)] ]

4. 高斯求积

前面几种算法一般要求小区间的宽度相同，即将区间等距划分，高斯法则不要求等距，其一般形式为：$$int^b_{a}f(x)dx approx sum^n_{i=1}C_if(x_i)$$ 其中(C_i)为点(x_i)对应的权重，点(x_i)是待定的高斯点。现在考虑([-1,1])区间：$$ int^1_{-1}f(x)dx approx sum^n_{i=1}C_i f(x_i)$$ 对(n=2)的情况，有：

[int^1_{-1}1dx=C_1+C_2=2 ]

[int^1_{-1}xdx=C_1x_1+C_2x_2=0 ]

[int^1_{-1}x^2dx=C_1x^2_1+C_2x^2_2=frac{2}{3} ]

[int^1_{-1}x^3dx=C_1x^3_1+C_2x^3_2=0 ]

这是一组非线性方程，存在无数个解，因此，额外提出对称性要求，即(x_1=-x_2)，则可求解得：$$ C_1=1,C_2=1,x_1=-frac{1}{sqrt{3}}=-0.57735027,x_2=frac{1}{sqrt{3}}=0.57735027 $$
即 $$int^1_{-1}f(x)dx approx f(-frac{1}{sqrt{3}})+ f(frac{1}{sqrt{3}})$$
同样地，可以推导出(n=3)的情况，当(n)越大时，计算的精度越高。下表给出了(n=2,3,4,5,6)时，对应的高斯点和权重：

上述推导过程是以([-1,1])区间为例，对积分区间([a,b])，可以通过如下线性变换，将积分区间变换为([-1,1])：$$x=frac{1}{2}[t(b-a)+a+b],dx = frac{1}{2}(b-a)dt$$

[int^b_{a}f(x)dx=int^1_{-1}f(frac{(b-a)t+a+b}{2})frac{(b-a)}{2}dt ]

5. 总结

本节课主要介绍函数积分的数值算法，最简单的有矩形法（左右点法）和中点法，在中点法的基础上，线性插值法提供了一种梯形求解算法，而采用非线性插值，则可以推导出Simpson 1/3 法和 3/8 法。这几种算法都要求对区间等距划分，最后介绍的高斯求积法则是确定高斯点和相应的权重来直接计算积分值。与数值微分类似，在数值积分算法中，也可以用Richardson外推加速算法将两个精度不高的方法组合构建出一个精度更高的数值积分算法，这里不做详细讨论，包括多重积分、多元积分等相关知识，这里暂不做深入研究，如果后面要使用，再进一步讨论。至此，暑期学校数值方法与计算课程的内容已经整理完毕。原本打算两周内整理好，后来因为各种事情耽搁了，以后尽量避免长线作战，趁热打铁方能加深印象。