数值计算笔记之数值积分(一)

2023-11-13

0、引言

一、数值积分的积分思想

1、中矩形公式

0、引言

在高数中，可以根据 $\int_{a}^{b}f(x)dx = F(b)-F(a)$ ，求得积分。但是如果 $f(x)$ 存在以下两种情况：

$f(x)$ 是离散函数
无法求得原函数 $F(x)$

那么对于这种情况，就可以利用数值积分。

所谓数值积分，指利用被积函数 $f(x)$ 在有限个点上的函数值来计算积分近似值的一种方法。

一、数值积分的积分思想

积分中值定理： $\int_{a}^{b}f(x)dx = (b-a)f(\xi ),\xi \in [a,b]$

一重积分可以理解为求面积，积分中值定理可以看成：在 $x\in [a,b]$ 上，函数 $f(x)$ 与 $x$ 轴的面积等于以 $(b-a)$ 为长， $\xi \in [a,b]$ 中一点函数值 $f(\xi )$ 为宽的矩形的面积。

1、中矩形公式

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

2、梯形公式

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

3、辛普森公式

$\int_{a}^{b}f(x)dx\approx (b-a)\cdot \frac{1}{6}[f(a)+4f(\frac{a+b}{2})+f(b)]$

一般情况下，可用： $f(\xi )\approx \sum_{k=0}^{n}w_{k}\cdot f(x_{k})$

$\Rightarrow \int_{a}^{b}f(x)dx = (b-a)\cdot f(\xi )\approx (b-a)\cdot \sum_{k=0}^{n}w_{k}\cdot f(x_{k}) = \sum_{k=0}^{n}A_{k}f(x_{k})$ ，其中 $A_{k}=(b-1)w_{k}$ 。

上式称为数值积分的基本公式， $x_{k}$ 称为求积节点， $A_{k}$ 称为求积系数。

二、求积公式的余项和代数精度

1、余项

$R\begin{bmatrix} f \end{bmatrix}=\int_{a}^{b}f(x)dx-\sum_{k=0}^{n}A_{k}\cdot f(x_{k})$

2、代数精度

若数值积分公式对任意的 $\leq m$ 次的代数多项式都准确成立，而对于 $x^{m+1}$ 却不能成立，则称该数值积分公式的代数精度为 $m$ 。

举个例子理解一下：求梯形公式的代数精度。

解： $\int_{a}^{b}f(x)dx\approx (b-a)\cdot \frac{f(a)+f(b)}{2}$

当 $f(x)=1$ ， $left = b-a = right$
当 $f(x)=x$ ， $left = \frac{1}{2}(b^{2}-a^{2})= right$
当 $F(x)=x^{2}$ ， $left =\frac{1}{3}(b^{3}-a^{3})$ ； $right = \frac{b-a}{2}(a^2+b^{2})$ $\therefore left \neq right$

$\therefore$ 梯形公式的代数精度为1。

定理：含有 n+1 个节点的插值型数值积分公式的代数精度至少为 n 。

三、插值型求积公式

已知

$x$	$x_{0}$	$x_{1}$	$\cdots$	$x_{n}$
$f(x)$	$f(x_{0})$	$f(x_{1})$	$\cdots$	$f(x_{n})$

则拉格朗日插值可表示为： $L_{n}(x)=\sum_{k=0}^{n}f(x_{k})\cdot l_{k}(x)$ ， $l_{k}(x) = \prod_{j=0;j\neq k}^{n}\frac{x-x_{j}}{x_{k}-x_{j}}$ 为基函数。

$\therefore \int_{a}^{b}f(x)dx \approx \int_{a}^{b}L_{n}(x)\cdot dx$

$= \int_{a}^{b}\sum_{k}^{n}f(x_{k})\cdot l_{k}(x)\cdot dx$

$= \sum f(x_{k})\cdot \int_{a}^{b}l_{k}(x)dx$

$= \sum A_{k}\cdot f(x_{k})$

即为插值型求积公式。其中 $A_{k}=\int_{a}^{b}l_{k}(x)dx$

四、牛顿--柯特斯公式 (N-C公式)

牛顿-柯斯特求积公式是插值型求积公式的特殊形式，在插值型求积公式中所取节点是等距时称为牛顿-柯斯特公式。

即，N-C公式为等距节点的插值型求积公式。

在[ $[a,b]$ 区间内设置等距的插值节点 $a=x0<x<\cdots <x_{n}=b$ .

设节点步长为 $h=\frac{b-a}{n}$

则节点 $x_{k}=a+k\cdot h ,(k=0,1,\cdots ,n)$

$\forall x\in [a,b]$ ，有 $x = a+t\cdot h$ $(0\leq t\leq n)$

则插值型求积系数：

$A_{k}=\int_{a}^{b}l_{k}(x)dx =\int_{a}^{b}\prod_{j=0;j\neq k}^{n}\frac{x-x_{j}}{x_{k}-x_{j}}dx$

$= \int_{0}^{n}\prod_{j=0;j\neq k}^{n} \frac{(a+th)-(a+jh)}{(a+kh)-(a+jh)}\cdot hdt$ 变量由 $x\rightarrow t$ ，所以变量范围由 $[a,b]\rightarrow [0,n]$ ； $dx = h\cdot dt$

$=\int_{0}^{n}\prod_{j=0;j\neq k}^{n}\frac{t-j}{k-j}\cdot hdt$

$=\int_{0}^{n}\frac{1}{k(k-1)\cdots (k-(k-1))(k-(k+1))\cdots (k-n)}\prod_{j=0;j\neq k}^{n}(t-j)\cdot hdt$

$=\int_{0}^{n}\frac{1}{k(k-1)\cdots 1(-1)(-2)\cdots (-(n-k))}\prod_{j=0;j\neq k}^{n}(t-j)\cdot hdt$ ，与 $t$ 无关的提到积分号之外

$=\frac{h}{k(k-1)\cdots 1(-1)(-2)\cdots (-(n-k))}\int_{0}^{n}\prod_{j=0;j\neq k}^{n}(t-j)\cdot dt$

$=\frac{\frac{b-a}{n}}{k!(n-k)!}(-1)^{n-k}\int_{0}^{n}\prod_{j=0;j\neq k}^{n}(t-j)\cdot dt$

$=(b-a)\frac{(-1)^{n-k}}{k!(n-k)!\cdot n}\int_{0}^{n}\prod_{j=0;j\neq k}^{n}(t-j)\cdot dt$

因为 $\frac{(-1)^{n-k}}{k!(n-k)!\cdot n}\int_{0}^{n}\prod_{j=0;j\neq k}^{n}(t-j)\cdot dt$ 只与 $k,n$ 有关，与 $a,b,f(x)$ ，步长 $h$ 均无关

而该表达式称之为柯斯特求积系数，记为

$C_{k}^{(n)}=\frac{(-1)^{n-k}}{k!(n-k)!\cdot n}\int_{0}^{n}\prod_{j=0;j\neq k}^{n}(t-j)\cdot dt$

$\therefore$ 求积系数 $A_{k} = (b-a)C_{k}^{(n)}$

N-C公式为： $\int_{a}^{b}f(x)dx = \sum_{k=0}^{n}A_{k}\cdot f(x_{k})=(b-a)\sum_{k=0}^{n}C_{k}^{(n)}\cdot f(x_{k})$

1、n = 1 (梯形公式)

$C_{k=0}^{(1)}=\frac{(-1)^{1-0}}{0!(1-0)!\cdot 1}\int_{0}^{1}(t-1)dt=\frac{1}{2}$

$C_{k=1}^{(1)}=\frac{(-1)^{1-1}}{1!(1-1)!\cdot 1}\int_{0}^{1}(t-0)dt=\frac{1}{2}$

$\therefore \int_{a}^{b}f(x)dx =(b-a)\cdot \frac{1}{2}\begin{bmatrix} f(a)+f(b) \end{bmatrix}$ 就变成了梯形公式

2、n = 2 (辛普森公式)

$C_{k=0}^{2} = \frac{(-1)^{2-0}}{0!(2-0)!\cdot 2}\int_{0}^{2}(t-1)(t-2)dt = \frac{1}{6}$ ， $C_{k=1}^{2} = \frac{4}{6}$ ， $C_{k=2}^{2} = \frac{1}{6}$

$\therefore \int_{a}^{b}f(x)dx = (b-a)\begin{bmatrix} \frac{1}{6}f(x_{0}+\frac{4}{6}f(x_{1})+\frac{1}{6}f(x_{2}) \end{bmatrix}$ $=\frac{(b-a)}{6} \begin{bmatrix} f(a)+4f(\frac{a+b}{2})+f(b) \end{bmatrix}$ 辛普森公式

3、n=4

$\int_{a}^{b}f(x)dx = \frac{b-a}{90}[7f(x_{0})+ 32f(x_{1})+12f(x_{2})+32f(x_{3})+7(x_{4})]$ ，其中 $x_{k} = a+k\cdot \frac{b-a}{4}$

五、复化求积法

当积分区间较大时，C-N求积公式的误差会很大。为减少误差，可以考虑增加节点，但用高次多项式插值导出的公式稳定性不好。所以可以考虑用分段插值函数来代替被积函数。

把 $[a,b]$ 分成n 个区间 $[x_{k-1},x_{k}] (k=1,2,\cdots ,n)$ ，然后给每个区间上用低阶N-C公式求积（设为 $I_{k}$ )，则

$\int_{a}^{b}f(x)dx \approx \sum_{k=1}^{n}I_{k}$

1、复化梯形公式

$T_{n}= \frac{h}{2}[f(x_{0})+f(x_{1})] + \frac{h}{2}[f(x_{1})+f(x_{2})] + \cdots + \frac{h}{2}[f(x_{n-1})+f(x_{n})]$

$= \frac{h}{2}[f(a)+ 2\sum_{k=1}^{n-1}f(x_{k}) + f(b)]$ ，其中 $h = \frac{b-a}{n}$

拓展：变步长梯形公式

定步长复化求积公式的一个明显缺点是：事先很难估计分划数n使结果达到预期精度。由于适当加密分点，精度会有所改善，为此采用自动加密分点的方法，并利用事后估计来控制加密次数，以判断是否达到预期精度，从而停止计算。

对所有已存在的子区间进行二分化，区间数由n变为2n
利用区间数为n时的积分值Tn以及新增的节点（即原来各子区间的中点）递推出区间数为2n时的积分值T2n
利用两次计算结果的差来估计误差，直到满足精度

公式如下：

$T_{2n} = \frac{1}{2}(T_{n}(h) + H_{n}(h))$ ，其中 $H_{n}(h) = h\sum_{k=0}^{n-1}f(a+(k+\frac{1}{2})h)$

2、复化辛普森公式（要求 n 为偶数）

$S_{n} = \frac{h}{3}[f(a) + 2\sum_{k=1}^{m-1}f(x_{2k}) + 4\sum_{k=1}^{m}f(x_{2k-1}) +f(b)]$ ,其中 $m = \frac{n}{2}$

六、代码实现

1、C++实现

先计算一个简单的，可以求原函数的。如 $f(x) = \int_{1}^{2}\frac{1}{x}$ ，标准答案为 $ln2 \approx 0.69314718$ 。程序运行结果为：

可以看出，复化公式的精度更高一些。

现在，如果要求原函数不可求的，如 $f(x) = \int_{1}^{2}\frac{log(3+x^{2})}{1+x^{2}}$ ,程序结果为：

C++代码:

main文件

#include<iostream>
#include"integral.h"
#include<iomanip>
using namespace std;

const double eps = 1e-8;

int main()
{
	double resultT = T(fun2,  1, 2);
	cout << "梯形公式计算结果" << "\t" << setprecision(8) << resultT << endl;

	double resultS = S(fun2, 1, 2);
	cout << "辛普森公式计算结果" << "\t" << setprecision(8) << resultS << endl;

	double resultTn = Tn(fun2, 1, 2, 100);
	cout << "复化梯形公式计算结果" << "\t" << setprecision(8) << resultTn << endl;

	double resultT2n = T2n(fun2, 1, 2, eps);
	cout << "变步长梯形公式计算结果" << "\t" << setprecision(8) << resultT2n << endl;

	cout << endl;
	double resultSn = Sn(fun2, 1, 2);
	cout << "复化辛普森公式计算结果" << "\t" << setprecision(8) << resultSn << endl;
}

.h头文件（积分实现部分）

#pragma once
#include<iostream>
using namespace std;
double fun1(double x) {
	double y = 1 / x;
	return y;
}
double fun2(double x) {
	double y = log(3 + pow(x,2)) / (1 + pow(x, 2));
	return y;
}


//梯形公式
double T(double(*f)(double), double a, double b) {
	double result = 0.0;
	result = 0.5 * (b - a) * (f(a) + f(b));
	return result;
}

//辛普森公式
double S(double(*f)(double), double a, double b) {
	double result = 0.0;
	result = 1 / 6.0 * (b - a) * ( f(a) + f(b) + 4 * f( (a + b) / 2.0) );
	return result;
}

//复化梯形公式
double Tn(double(*f)(double), double a, double b, int n) {
	double result = 0.0;
	double h = (b - a) / n;
	result = h / 2.0 * (f(a) + f(b));
	for (int k = 1; k < n; k++) {
		double x = a + k * h;
		result += h * f(x);
	}
	return result;
}

//变步长梯形积分
double T2n(double(*f)(double), double a, double b, double error) {
	int n = 1;   
	double h = (b - a) / n;
	double Tn = (f(a) + f(b)) * h / 2.0;  //计算 n=1 时的积分，粗略值
	double T2n;     //步长变为一半后的积分值
	bool key = false;
	do {    //至少循环一次
		double Hn = 0.0;
		for (int k = 0; k < n; k++) {
			double x = a + (k + 0.5) * h;
			Hn += h * f(x);
		}
		T2n = (Tn + Hn) / 2.0;
		//判断误差
		if (fabs(Tn - T2n) < error) {
			key = true;
		}
		else {
			Tn = T2n;
			n *= 2;
			h /= 2;
		}

	} while (!key);
	return T2n;
}

//复化辛普森公式
double Sn(double(*f)(double), double a, double b) {
	int n;
	cout << "正在使用复化辛普森公式，请输入分段个数n, 要求为偶数。 n =";
	cin >> n;
	/*cout << endl;*/
	while (n % 2) {
		cout << "输入n要求为偶数，请重新输入, n = ";
		cin >> n;
		cout << endl;
	}
	double result = 0.0;
	double h = (b - a) / n;
	result = h / 3.0 * (f(a) + f(b));
	for (int k = 1; k < n; k++) {
		double x = a + k * h;
		if (k % 2 == 0) {
			result += h / 3.0 * 2 * f(x);
		}
		else {
			result += h / 3.0 * 4 * f(x);
		}
	}
	return result;
}

2、matlab实现

matlab中自带了函数，例如， $quadtx$ 函数、 $quadgui$ 函数。

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)

数值计算

c