文章目录

前言

1 数值微分

2 数值积分

小结

前言

今天我们要说的就是数值微积分，赶紧看看他和高等数学中的微积分有什么区别吧。本文是科学计算与MATLAB语言专题六第一小节的学习笔记，如果大家有时间的话，可以去听听课，没有的话，可以看看我的笔记，还是很不错的。

1 数值微分

(1)数值差分与差商

微积分中，任意函数f ( x ) f(x)f(x)在x 0 x_0x0​点的导数是通过极限定义的：

f ′ ( x 0 ) = lim ⁡ t → 0 f ( x 0 + h ) − f ( x 0 ) h f'(x_0)=\lim\limits_{t \to 0}\frac{f(x_0+h)-f(x_0)}{h}f′(x0​)=t→0lim​hf(x0​+h)−f(x0​)​

f ′ ( x 0 ) = lim ⁡ t → 0 f ( x 0 ) − f ( x 0 − h ) h f'(x_0)=\lim\limits_{t \to 0}\frac{f(x_0)-f(x_0-h)}{h}f′(x0​)=t→0lim​hf(x0​)−f(x0​−h)​

f ′ ( x 0 ) = lim ⁡ t → 0 f ( x 0 + h / 2 ) − f ( x 0 − h / 2 ) h f'(x_0)=\lim\limits_{t \to 0}\frac{f(x_0+h/2)-f(x_0-h/2)}{h}f′(x0​)=t→0lim​hf(x0​+h/2)−f(x0​−h/2)​

如果去掉极限定义中h趋向于0的极限过程，得到函数在x 0 x_0x0​点处以h(h>0)为步长的向前差分、向后差分和中心差分公式：

向前差分：Δ r ( x 0 ) = f ( x 0 + h ) − f ( x 0 ) \Delta r(x_0)=f(x_0+h)-f(x_0)Δr(x0​)=f(x0​+h)−f(x0​)

Δ \DeltaΔ原来她叫Delta

向后差分：∇ P ( x 0 ) = f ( x ) − f ( x 0 − h ) \nabla P(x_0)=f(x)-f(x_0-h)∇P(x0​)=f(x)−f(x0​−h)

∇ \nabla∇她叫nabla

中心差分：δ F ( x 0 ) = f ( x 0 + h / 2 ) − f ( x 0 − h / 2 ) \delta F(x_0)=f(x_0+h/2)-f(x_0-h/2)δF(x0​)=f(x0​+h/2)−f(x0​−h/2)

δ \deltaδ她叫delta

当步长h充分小时，得到函数在x 0 x_0x0​点处以h(h>0)为步长的向前差商、向后差商和中心差商公式：

向前差分：f ′ ( x 0 ) ≈ f ( x 0 + h ) − f ( x 0 ) f'(x_0)\approx f(x_0+h)-f(x_0)f′(x0​)≈f(x0​+h)−f(x0​)

向后差分：f ′ ( x 0 ) ≈ f ( x ) − f ( x 0 − h ) f'(x_0)\approx f(x)-f(x_0-h)f′(x0​)≈f(x)−f(x0​−h)

中心差分：f ′ ( x 0 ) ≈ f ( x 0 + h / 2 ) − f ( x 0 − h / 2 ) f'(x_0)\approx f(x_0+h/2)-f(x_0-h/2)f′(x0​)≈f(x0​+h/2)−f(x0​−h/2)

函数f(x)在点x0的微分接近于函数在该点的差分，而f在点x的导数接近于函数在该点的差商。

(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，按行计算差分。

注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样，对于矩阵来说，差分后的矩阵比原矩阵少了一行或一列。

另外，计算差分之后，可以用f(x)在某点处的差商作为其导数的近似值。

例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.0456

2 数值积分

2.数值积分

(1)数值积分基本原理

在高等数学中，计算定积分依靠微积分基本定理，只要找到被积函数f(x)

的原函数大(x)，则可用牛顿一莱布尼兹(Newton-Leibniz)公式：

∫ a b f ( x ) d x = F ( b ) − F ( a ) \int _a^bf(x)dx=F(b)-F(a)∫ab​f(x)dx=F(b)−F(a)在有些情况下，应用牛顿一莱布尼兹公式有困难，例如，当被积函数的原函数无法用初等函数表示，或被积函数是用离散的表格形式给出的。这时就需要用数值解法来求定积分的近似值。

求定积分的数值方法多种多样，如梯形法、辛普森(Simpson)法、高斯求积公式等。它们的基本思想都是将积分区间[a，b]分成n个子区间[ x i ， x i + 1 ] [x_i，x_{i+1}][xi​，xi+1​]，i=1，2，…，n，其中x i = a x_i=axi​=a，x i + 1 = b x_{i+1}=bxi+1​=b，这样求定积分问题就分解为下面的求和问题：

S = ∫ a b f ( x ) d x = ∑ i = 1 n ∫ x i x i + 1 f ( x ) d x S= \int_a^b f(x)dx=\sum\limits _{i=1}^n\int_{x_i}^{x_{i+1}} f(x)dxS=∫ab​f(x)dx=i=1∑n​∫xi​xi+1​​f(x)dx在每一个小的子区间上定积分的值可以近似求得，从而避免了牛顿一莱布尼兹公式需要寻求原函数的困难。

(2)数值积分的实现

基于自适应辛普森方法

[I，n]=quad(filename，a，b，tol，trace)

基于自适应Gauss-Lobatto方法

[I，n]=quad1(filename，a，b，tol，trace)

其中，filename是被积函数名；

a和b分别是定积分的下限和上限，积分限[a，b]必须是有限的，不能为无穷大(Inf)；

tol用来控制积分精度，默认时取tol=1 0 − 6 10^{-6}10−6；

trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0；返回参数I即定积分的值，n为被积函数的调用次数。

例2 分别用quad函数和quadl函数求定积分∫ 0 1 1 1 + x 2 d x \int_0^1\frac{1}{1+x^2}dx∫01​1+x21​dx的近似值，并在相同的积分精度下，比较被积函数的调用次数。

format long%有效数字保留16位。

f=@(x) 4./(1+x.^2);

[I,n]=quad(f,0,1,1e-8)%基于自适应辛普森方法求近似值。

[I,n]=quadl(f,0,1,1e-8)%基于自适应Gauss-Lobatto方法求近似值。

(atan(1)-atan(0))*4%求理论值。

format short%format short：恢复默认格式，小数点后保留4位。

I =

3.141592653733437

n =

61

I =

3.141592653589806

n =

48

ans =

3.141592653589793

基于全局自适应积分方法

I=integral(filename，a，b)

其中，I是计算得到的积分；

filename是被积函数；

a和b分别是定积分的下限和上限，积分限可以为无穷大。

例3 求定积分I = ∫ 1 e f = 1 ( x ( 1 − ln ⁡ 2 x ) ) I=\int_1^ef=\frac{1}{(x\sqrt{(1-\ln^2x))}}I=∫1e​f=(x(1−ln2x))​1​。

被积函数文件fe.m：

function f=fe(x)

f=1./(x.*sqrt(1-log(x).^2));

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

基于自适应高斯-克朗罗德方法

[I，err]=quadgk(filename，a，b)

其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大(-Inf或Inf)，也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。

例4 求定积分∫ 2 π + ∞ f = 1 x 2 s i n 1 x d x \int _\frac{2}{\pi}^{+\infty }f=\frac{1}{x^2}sin\frac{1}{x}dx∫π2​+∞​f=x21​sinx1​dx。

被积函数文件fsx.m：

function f=fsx(x)

f=sin(1./x)./x.^2;

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

I =

1.0000

基于梯形积分法

已知( x i ， y i ) ( i = 1 ， 2 ， … ， n ) (x_i，y_i)(i=1，2，…，n)(xi​，yi​)(i=1，2，…，n)，且a = x 1 < x 2 < … < x n = b a=x_1a=x1​

I=trapz(x，y)

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

trapz函数采用梯形积分法则，积分的近似值为：

I = ∑ i = 1 n − 1 h i y i + 1 + y i 2 I=\sum\limits_{i=1}^{n-1}h_i\frac{y_{i+1}+y_i}{2}I=i=1∑n−1​hi​2yi+1​+yi​​其中，h i = x i + 1 − x i h_i=x_{i+1}-x_ihi​=xi+1​−xi​。

可用以下语句实现：

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)

(3)多重定积分的数值求解

求二重积分的数值解：∫ c d ∫ a b f ( x ， y ) d x d y \int _c^d \int_a^b f(x，y)dxdy∫cd​∫ab​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)

求三重积分的数值解：∫ e f ∫ c d ∫ a b f ( x ， y ， z ) d x d y d z \int_e^f \int_c^d \int_a^b f(x，y，z)dxdydz∫ef​∫cd​∫ab​f(x，y，z)dxdydz

I=integra13(filename，a，b，c，d，e，f)

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

例6 分别求二重积分和三重积分。

∫ − 1 1 ∫ − 2 2 e − x 2 / 2 s i n ( x 2 + y ) d x d y \int_{-1}^1\int_{-2}^2e^{-x^2/2 }sin(x^2+y)dxdy∫−11​∫−22​e−x2/2sin(x2+y)dxdy

∫ 0 1 ∫ 0 π ∫ 0 π 4 x z e − z 2 y − x 2 d x d y d z \int_0^1\int_0^\pi\int_0^\pi 4xze^{-z^2y-x^2}dxdydz∫01​∫0π​∫0π​4xze−z2y−x2dxdydz

f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);

I1=quad2d(f1,-2,2,-1,1)

f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x);

I2=integral3(f2,0,pi,0,pi,0,1)

小结

大家学会了吗？最后欢迎大家

点赞，收藏⭐，转发，

如有问题、建议，请您在评论区留言哦。