数值分析(Matlab): 数值微分

Amor灬Niss

于 2023-03-08 21:25:50 发布

阅读量1k

点赞数

分类专栏：数值分析文章标签： matlab

本文链接：https://blog.csdn.net/qq_37896521/article/details/129411127

版权

数值分析专栏收录该内容

2 篇文章 0 订阅

订阅专栏

文章介绍了差商法用于求解一阶和二阶导数的过程，特别讨论了自变量为矢量的函数。一阶导数的一般函数和矢量函数的差商公式被详细阐述，包括向前、向后和中心差商。二阶导数在矢量函数中对应于Jacobian矩阵，文章给出了中心差分的数值算法来计算它。

摘要由CSDN通过智能技术生成

前言

(无)

1. 差商法求一阶导数

1.1 一般函数

由导数定义:
${{{\rm{d}}f} \over {{\rm{d}}x}} = \mathop {\lim }\limits_{h \to 0} {{f\left( {x + h} \right) - f\left( x \right)} \over h}$
则可以直接使用该式进行求导计算，称之为差商法求导。同时用Taylor公式展开，可以得到其余项公式，即数值精度：

向前差商公式
${{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( {x + h} \right) - f\left( x \right)} \over h} + o\left( h \right)$
向后差商公式
${{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( x \right) - f\left( {x - h} \right)} \over h} + o\left( h \right)$
中心差商公式
${{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( {x + h} \right) - f\left( {x - h} \right)} \over {2h}} + o\left( {{h^2}} \right)$

1.2 自变量为矢量的函数

与1.1节中的公式一样，只是需要对自变量的分量逐个差商。

% 假设   U = x * y^2 * z^3
% 则     F = dU/dx = [y^2*z^3; x*z^3*2y; xy^2*3z^2];
%        J = [  0,         2y*z^3,    3y^2*z^2;
%               2*z^3*y,   2x*z^3,    6x*y*z^2;
%               3*y^2*z^2, 6x*y*z^2,  6xy^2z];
x = [1.0; 2.0; 3.0];
h = 1e-5;
F = zeros(3,1);
e = cell(3);
e{1} = [1;0;0];
e{2} = [0;1;0];
e{3} = [0;0;1];
for i=1:3
    F(i) = 1/(2*h) * (...
               Ep(x + h*e{i})...
             - Ep(x - h*e{i})...
        );
end

F_ans = [ x(2)^2 * x(3)^3;...
    2* x(1) * x(3)^3 * x(2);...
    3* x(1) * x(2)^2 * x(3)^2];
Dif_F = [F, F_ans, F-F_ans];

%%
function U = Ep(x)
U = x(1) * x(2)^2 * x(3)^3;
end

结果为：
在这里插入图片描述

2. 二阶导数

2.1 自变量为矢量

自变量为矢量时，二阶导数通常称为Jacobian，是一个矩阵。其中心差分可以表示为：
在这里插入图片描述
其中：
${\bf{q}} = \sum\limits_i {{q_i}{{\rm{e}}_i}}$
对应的数值算法为：

% 假设   U = x * y^2 * z^3
% 则     F = dU/dx = [y^2*z^3; x*z^3*2y; xy^2*3z^2];
%        J = [  0,         2y*z^3,    3y^2*z^2;
%               2*z^3*y,   2x*z^3,    6x*y*z^2;
%               3*y^2*z^2, 6x*y*z^2,  6xy^2z];
x = [1.0; 2.0; 3.0];
h = 1e-5;
e = cell(3);
e{1} = [1;0;0];
e{2} = [0;1;0];
e{3} = [0;0;1];
%%
a = zeros(3,3);
x = [1.0; 2.0; 3.0];
for i = 1:3
    for j = 1:3
        a(i,j) = 0.25/(h^2) * (...
              Ep(x + h*e{i} + h*e{j})...
            - Ep(x + h*e{i} - h*e{j})...
            - Ep(x - h*e{i} + h*e{j})...
            + Ep(x - h*e{i} - h*e{j})...
        );
    end
end
J = [  0,               2*x(2)*x(3)^3,      3*x(2)^2*x(3)^2;
       2*x(3)^3*x(2),   2*x(1)*x(3)^3,      6*x(1)*x(2)*x(3)^2;
       3*x(2)^2*x(3)^2, 6*x(1)*x(2)*x(3)^2, 6*(1)*x(2)^2*x(3)];
J = reshape(J, 9, 1);
a = reshape(a, 9, 1);
Dif_J = [J, a, J-a];

%%
function U = Ep(x)
U = x(1) * x(2)^2 * x(3)^3;
end