matlab 数值微分和积分

数值微分和积分

数值微分

（1）数值微分和差商

向前差商：
$$f^{{}^{\prime }}(x_0)=\frac{f(x_0+h)-f(x_0)}{h}$$
向后差商：
$$f^{{}^{\prime }}(x_0)=\frac{f(x_0)-f(x_0-h)}{h}$$
中心差商：
$$f^{{}^{\prime }}(x_0)=\frac{f(x_0+\frac{h}{2})-f(x_0-\frac{h}{2})}{h}$$

（2）数值微分的实现

MATLAB提供了求向前差分的函数 diff，其调用格式有三种

• dx = diff(x)：计算向量x的一阶向前差分，dx(i) = x(i+1)-x(i),i=1,2,……，n-1。
• dx = diff(x,n)：计算向量 x 的 n阶向前差分。例如diff(x,2) = diff(diff(x))。
• dx = diff(A,n,dim): 计算矩阵A的 n 阶差分，dim = 1(默认状态），按列计算差分；dim = 2，按行计算差分。

例1：设 f(x) = sinx，在[0,2π]范围内随机采样，计算 f '(x) 的近似值，并与理论值 f '(x)=cosx进行比较。

x = [0,sort(2*pi*rand(1,5000),2*pi)]
y = sin(x);
f1 = diff(y)./diff(x)
f2 = cos(x(1:end-1))
plot(x(1:end-1),f1,x(1:end-1),f2);
d = norm(f1-f2)			% d= 0.0433，求范数，相接近

数值积分

(1)数值积分基本原理

牛顿-莱布尼兹（Newton-Leibniz）公式：
$${\int }_a^{b}f(x)dx=F(b)-F(a)$$
将积分区间[a,b]分成 n 个子区间 [xi, x_i+1]，i = 1,2,……,n，其中x1= a, x_(n+1) = b 这样求定积分问题就分解为下面的求和问题
$$S={\int }_a^{b}f(x)dx=\sum {\int }_{x_i}^{x_{i+1}}f(x)dx$$

(2)数值积分的实现

• 基于自适应辛普森方法

  [l , n] = quad(filename, a, b, tol, trace)

• 基于自适应Gauss-Lobatto 方法

  [l , n] = quadl (filename, a, b, tol, trace)

  其中，

  1、filename是被积函数名；

  2、a和b分别是定积分的下限和上限，积分限[a, b] 必须是有限的，不能为无穷大( lnf );

  3、tol 用来控制积分精度，默认时取 tol = 10^-6;

  4、 trace 控制是是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时 trace = 0；

  5、返回参数 l 即定积分的值, n为被积函数的调用次数。

例2 分别用quad 和 quadl 函数求定积分的近似值，并在相同积分精度下，计较被积函数调用的次数。
$${\int }_0^1\frac{4}{1+x^2}dx$$

format long
f = @(x) 4./(1+x.^2)
[l,n] = quad(f,0,1,1e-8)
% l =  3.141592653733437
% n = 61

[l,n] = quadl(f,0,1,1e-8)
%l = 3.141592653589806
%n = 48

• 基于全局自适应积分方法

  l = integral(filename, a, b)

其中

1、l 是计算得到的积分；

2、filename 是被积函数；

3、a和b是积分上下限，积分限可以为无穷大

例3 求定积分
$$I={\int }_1^e\frac{1}{x\sqrt{1-l{n}^{2}x}}dx$$

%被积函数文件 fe.m
function f = fe(x)
f = 1./(x.*sqrt(1-log(x).^2));


>>I = integral(@fe,1,exp(1))
I = 1.5708

• 基于自适应高斯-克朗罗德方法**（计算振荡函数的定积分）**

  [l, err] = quadgk(filename,a,b)

  1、err 返回近似误差范围

  2、其他参数的含义和用法与 quad 函数相同。

  3、积分上下限可以是无穷大（-lnf 或 lnf），也可以是复数，如果积分上下限是复数，则quadgk函数在复平面上求积分

例4、求定积分
$${\int }_{\frac{2}{\pi }}^{+\infty }\frac{1}{x^2}sin\frac{1}{x}dx$$

% 被积函数文件 fsx.m
function f = fsx(x)
f = sin(1./x)./x.^2

>> I = quadgk(@fsx,2/pi,+Inf)

I =

    1.0000

• 基于梯形积分法（当不知道/被积函数，只知道部分积分点时）

  已知（xi, yi) (i = 1,2,……,n)，且a = x1<x2<x3<……<xn = b

  求 I = ${\int }_a^{b}f(x)dx$ 的近似值。

  I = trapz(x,y)

  其中，向量x, y 定义函数关系 y = f(x)。

  trapz 函数采用梯形积分法则，积分近似值为
  $$I=\sum _{i=1}^{n-1}hi\frac{y_{i+1}+y_i}{2}$$
  其中，hi = $x_{i+1}$-$x_i$

可以用以下语句实现

sum(diff(x).*(y(1:end-1)+y(2:end))/2)

例5 设x = 1:6, y = [6,8,11,7,5,2]，用trapz函数计算定积分。

x = 1:6;
y = [6,8,11,7,5,2];
plot(x,y,'-ko');
grid on
axis([1,6,0,11])
I1 = trapz(x,y)
I2 = sum(diff(x).*(y(1:end-1)+y(2:end))/2)
% I1 = I2 = 35

• 多重积分的数值求解

  1、求二重积分的数值解 ${\int }_c^d{\int }_a^{b}f(x,y)dxdy$

  I = integral2(filename,a,b,c,d)

  I = quad2d(filename,a,b,c,d)

  I = dblquad(filename,a,b,c,d,tol)

• 求三重积分的数值解 ${\int }_e^f{\int }_c^d{\int }_a^{b}f(x,y)dxdy$

  I = integral3(filename,a,b,c,d,e,f)

  I = triplequad(filename,a,b,c,d,e,f,tol)