【Matlab】数值微积分与方程求解

数值微积分与方程求解 数值微分

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

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

x=[0,sort(2*pi*rand(1,5000)),2*pi]; % 确保初始为0，末尾为2Π而已
y=sin(x);
f1=diff(y)./diff(x); % 差分模拟实现求导
f2=cos(x(1:end-1));
subplot(121)
plot(x(1:end-1),f2, 'r'); % 数值解
subplot(122)
plot(x(1:end-1),f1, 'g'); % 解析解
d=norm(f1-f2) % 欧氏距离，即平方和开根

数值积分

数值积分基本原理：牛顿-莱布尼兹公式$\int_a^b f(x),dx =F(b)-F(a) $

数值积分的实现：

• 基于自适应辛普森方法

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

• 基于自适应Gauss-Lobatto方法

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

  参数解释：

  filename是被积函数名或文件名或句柄；a和b分别是定积分的下限和上限，积分限[a, b]必须是有限的，不能为无穷大( lnf ) ；tol用来控制积分精度，默认时取tol=10e-6；trace控制是否展现积分过程，若取非0则展现积分过程，取0则不展现，默认时取trace=0；返回参数I即定积分的值，n为被积函数的调用次数。

求 I = ∫ 0 1 4 1 + x 2   d x I=\int_0^1 \frac{4}{1+x^2}\,dx I=∫01​1+x24​dx

format long;
f = @(x) 4./(1+x.^2);
[I, n]=quad(f, 0, 1, 1e-8) % quad函数
[I, n]=quadl(f, 0, 1, 1e-8) % quadl函数
(atan(1) - atan(0))*4 % 直接求解，解析解

结果如下：

I =

   3.141592653733437

n =

    61 % 调用quad时被积函数执行次数

I =

   3.141592653589806

n =

    48 % 调用quadl时被积函数执行次数，因此此函数执行效率高

ans =

   3.141592653589793

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

  I = integral(filename, a, b)

  参数解释：

  l是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。

求定积分 I = ∫ 1 e 1 x 1 − l n 2 x   d x I = \int_1^e \frac{1}{x\sqrt{1-ln^2x}}\,dx I=∫1e​x1−ln2x ​1​dx

I = integral(@(x)1./(x.*sqrt(1-log(x).^2)), 1, exp(1)); % 1.570796326795570

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

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

• 基于梯形积分法

  应用场景：已知样本点，求样本点构成函数的积分

  调用格式：I = trapz(x, y)，向量x、y定义函数关系 y=f(x)

  实现原理：

  ​ 已知 ( x i , y i )   ( i = 1 , 2 , . . . , n ) (x_i, y_i) \space (i=1,2,...,n) (xi​,yi​) (i=1,2,...,n)，且 a = x 1 < x 2 < . . . < x n = b a=x_1