前言
(无)
1. 差商法求一阶导数
1.1 一般函数
由导数定义:
d
f
d
x
=
lim
h
→
0
f
(
x
+
h
)
−
f
(
x
)
h
{{{\rm{d}}f} \over {{\rm{d}}x}} = \mathop {\lim }\limits_{h \to 0} {{f\left( {x + h} \right) - f\left( x \right)} \over h}
dxdf=h→0limhf(x+h)−f(x)
则可以直接使用该式进行求导计算,称之为差商法求导。同时用Taylor公式展开,可以得到其余项公式,即数值精度:
- 向前差商公式
d f d x = f ( x + h ) − f ( x ) h + o ( h ) {{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( {x + h} \right) - f\left( x \right)} \over h} + o\left( h \right) dxdf=hf(x+h)−f(x)+o(h) - 向后差商公式
d f d x = f ( x ) − f ( x − h ) h + o ( h ) {{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( x \right) - f\left( {x - h} \right)} \over h} + o\left( h \right) dxdf=hf(x)−f(x−h)+o(h) - 中心差商公式
d f d x = f ( x + h ) − f ( x − h ) 2 h + o ( h 2 ) {{{\rm{d}}f} \over {{\rm{d}}x}} = {{f\left( {x + h} \right) - f\left( {x - h} \right)} \over {2h}} + o\left( {{h^2}} \right) dxdf=2hf(x+h)−f(x−h)+o(h2)
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,是一个矩阵。其中心差分可以表示为:
其中:
q
=
∑
i
q
i
e
i
{\bf{q}} = \sum\limits_i {{q_i}{{\rm{e}}_i}}
q=i∑qiei
对应的数值算法为:
% 假设 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