数值积分及复变函数

中点公式、 Newton-Cotes(牛顿-科特斯)公式、 Gauss公式、三角形上的求积公式，以及MATLAB中提供的求积函数。

1 数值积分

本节主要讨论 $I(f)={\int }_a^{b}f(x)dx$ 的数值方法，$f$ 是区间 $[a,b]$ 上的任意连续函数。

1.1 中点公式

数值积分的基本问题就是针对某些函数类，选择合适的求积节点和求积函数，使得求积公式尽可能快地逼近 $I(f)$。

用区间中点的函数值近似代替该区间上的函数值，即 $I_{mp}^c={\sum }_{k=1}^M(x_k-x_{k-1})f(\frac{x_{k-1}+x_k}{2})$。该求积公式称为复合中点公式，该求积公式关于步长具有二阶精度。

常见的中点公式是复合中点公式的特殊情况，即当 $M=1$ 时的情况为 $I_{mp}=(b-a)f(\frac{a+b}{2})$。中点公式的MATLAB实现程序如下：

function y=mpint(a,b,M,f)
h=(b-a)/M;
x=linspace(a+h/2,b-h/2,M);
y=sum(feval(f,x)*h);

1.2 Newton-Cotes公式

由于多项式是最简单的函数类，构造求积公式的一个最基本构思是利用拉格朗日插值多项式。

复化梯形公式 $T_n=\frac{h}{2}[f(a)+f(b)+2{\sum }_{i=1}^{n-1}f(a+ih)]$，复化辛普森公式 $S_n=\frac{h}{6}{\sum }_{i=0}^{n-1}[f(x_i)+4f(x_{i+\frac{1}{2}})+f(x_{i+1})]$，其中，$x_{i+\frac{1}{2}}$ 表示区间 $[x_i,x_{i+1}]$ 的中点。

复化梯形公式和复化辛普森公式的MATLAB程序如下：

% 复化梯形公式
function y=traint(a,b,n,f)
h=(b-a)/n;
x=linspace(a,b,n+1);
y1=h*feval(f,x);
y1(1)=y1(1)/2;
y1(n+1)=y1(n+1)/2;
y=sum(y1);

% 复化辛普森公式
function y=sraint(a,b,n,f)
h=(b-a)/n;
x=linspace(a,b,2*n+1);
y1=feval(f,x);
y1(2:2:2*n)=4*y1(2:2:2*n);
y1(3:2:2*n-1)=2*y1(3:2:2*n-1);
y=h/6*sum(y1);

1.3 Guass求积公式

取 $n$ 次正交多项式的零点作为节点，则相应的求积公式具有 $2n-1$ 次代数精度。这种公式称为 $Gauss$ 型求积公式，并称其节点为 $Gauss$ 节点。

对于一般的区间 $[a,b]$，只要作一个简单的变换即可将 $[-1,1]$ 上的 $Gauss$ 积分点转换为区间 $[a,b]$ 上的 $Gauss$ 积分点。

$Gauss$ 求积公式的MATLAB程序如下：

function y = gauss_int(a,b,n,f)
switch n
    case 1
        x = 0;
        w = 2;
    case 2
        x = [0.5773502692;-0.577350269];
        w = [1;1];
    case 3
        x = [0.7745966692;-0.7745966692;0];
        w = [0.5555555556;0.5555555556;0.8888888889];
    case 4
        x = [0.8611363116;-0.8611363116;0.3399810436;-0.3399810436;];
        w = [0.3478548451;0.3478548451;0.6521451549;0.6521451549;];
    case 5
        x = [0.9061798459;-0.9061798459;0.5384693101;-0.5384693101;0];
        w = [0.2369268851;0.2369268851;0.4786286705;...
            0.4786286705;0.5688888889];
    case 6
        x = [0.9324695142;-0.9324695142;0.6612093865;...
            -0.6612093865;0.2386191861;-0.2386191861];
        w = [0.1713244924;0.1713244924;0.3607615730;...
            0.3607615730;0.4679139346;0.4679139346];
    case 7
        x = [0.9491079123;-0.9491079123;0.7415311856;...
            -0.7415311856;0.4058451514;-0.4058451514;0];
        w = [0.1294849662;0.1294849662;0.2797053915;...
            0.2797053915;0.3818300505;0.3818300505;0.417951837];
    case 8
        x = [0.9602898565;-0.9602898565;0.7966664774;-0.7966664774;...
            0.5255324099;-0.5255324099;0.1834346425;-0.1834346425];
        w = [0.1012285363;0.1012285363;0.2223810345;0.2223810345;...
            0.3137066459;0.3137066459;0.3626837834;0.3626837834];
    otherwise
        disp('n太大了');
        y = Inf;
        return;
end
x = (b - a) / 2 * x + (a + b) / 2;
y = sum(w .* feval(f,x)) * (b - a) / 2;

1.4 三角形上的求积公式

在实际计算当中，经常会用到三角形上求数值积分。如果三角形上的积分公式对三角形上的 $n$ 次多项式精确成立，就称该求积公式具有 $n$ 次代数精度。

三角形上求积公式MATLAB程序如下所示：

function y = triangl_int(a1,a2,a3,n,f)
s = abs(det([a1 1;a2 1;a3 1])) / 2;
switch n
	case 1
		y = s * feval(f,(a1+a2+a3)/3);
	case 2
		y = s / 3 * (feval(f,(a1+a2)/2) + feval(f,(a2+a3)/2)...
			+ feval(f,(a3+a1)/2));
	case 3
		y = s / 60 * 3*(feval(f,a1) + feval(f,a2) + feval(f,a3)) +...
			8*(feval(f,(a1+a2)/2) + feval(f,(a2+a3)/2) + ...
			feval(f,(a3+a1)/2)) + 27*feval(f,(a1+a2+a3)/3);
	otherwise
		disp('n 太大了！')
		y = Inf;
		return;
end

1.5 MATLAB提供的求积函数

1.5.1 int积分函数

int积分函数可用于求不定积分和定积分，其调用格式如下所示：

int(S)
int(S,v)
int(S,a,b)
int(S,v,a,b)	% 符号函数 S 定义的函数关于变量 v 从 a 到 b 的定积分

1.5.2 trapz函数

trapz函数是梯形数值积分公式，其调用格式如下所示：

Q=trapz(Y)
Q=trapz(X,Y)
Q=trapz(...,dim)

1.5.3 quad函数

quad函数可以对 M文件建立的函数作数值积分，采用的是自适应辛普森方法，属于低阶方法，精度较低。其调用格式如下所示：

q=quad(fun,a,b)
q=quad(fun,a,b,tol)
q=quad(fun,a,b,tol,trace)
[q,fcnt]=quad(...)

1.5.4 quadl函数

quad函数可以对 M文件建立的函数作数值积分，采用的是自适应Lobatto方法，属于高阶方法，精度较高。其调用格式如下所示：

q=quadl(fun,a,b)
q=quadl(fun,a,b,tol)
q=quadl(fun,a,b,tol,trace)
[q,fcnt]=quadl(...)

1.5.5 其他函数

除了上述求积函数之外，MATLAB还提供了自适应Gauss-Kronrod求积函数quadgk、矢量值函数的积分函数quadv、平面区域上的积分函数quad2d、矩形区域上的积分函数dblquad、三维矩阵上的积分函数triplequad等。

2 复变函数

2.1 复变函数的极限求导和积分

limit函数可用于求解复变函数的极限，也可用于求解实变函数的极限。其调用格式如下所示：

limit(F,x,a)	% 符号函数 F 当 x 趋向于 a 时的极限
limit(F,a)
limit(F)
limit(F,x,a,'right')	% 右极限
limit(F,x,a,'left')	% 左极限

diff函数可用于求解复变函数的导数，也可用于求解实变函数的导数。其调用格式如下所示：

Y=diff(X)
Y=diff(X,n)
Y=diff(X,n,dim)	% 符号表达式 X 在制定维度的 n 阶导数

复变函数的积分和实变函数的积分是一样的，区别在于复变函数的积分中可能会出现复数。

2.2 复变函数的Taylor展开

MATLAB中函数的Taylor级数展开由函数taylor实现，其调用格式如下所示：

taylor(f)
taylor(f,x)
taylor(f,x,a)	% 返回函数 f 在 a 点附近的五次幂多项式近似

以求解函数 $f(x)=si{n}^{2}x$ 在 $x=0$ 的Taylor展开级数为例。其MATLAB程序如下所示：

clc
syms x
y = sin(x)^2+sin(x)+cos(2*x);
taylor(y,x,0)

输出结果如下所示：

ans = (sym)

x^5/120+x^4/3-x^3/6-x^2+x+1

2.3 复变函数图像

绘制复变函数的图像主要使用命令plot，关于plot的具体细节见 绘图与Simulink简介。

以复变函数 $f(t)=e^{it}+e^{-it}$ 图形为例。其MATLAB程序如下：

clc
t=-5:0.1:5;
plot(exp(i*t)+exp(-i*t))

其结果如下所示：

2.4 留数

在MATLAB中，留数的计算是由residue来实现的，其调用格式如下所示：

[r,p,k]=residue(b,a)	% 返回由两个多项式 b(s) 和 a(s) 的商的留数矢量、极点矢量，以及两个多项式比值的部分分式展开的直接项，矢量 a 和 b 构成多项式 a(s)、b(s) 的系数
[b,a]=residue(r,p,k)	% 将部分分式展开成两多项式的商的形式，它们的系数放在矢量 b 和 a 中

参考文献

MATLAB从入门到实践（第2版）