文章目录一、问题提出二、一阶导数和二阶导数的二阶中心差分格式三、一阶导数和二阶导数的四阶中心差分格式四、MATLAB代码 一、问题提出 ,且横坐标等间距(即 且 时,恒有 ),如何计算/估计一阶导数 及二阶导数 ? 二、一阶导数和二阶导数的二阶中心差分格式![在这里插入图片描述 python二阶中心差法 二阶中心差分推导_一阶导数的四阶中心差分格式_07](https://s2.51cto.com/images/blog/202310/26162531_653a227bb182820062.png?x-oss-process=image/watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_30,g_se,x_10,y_10,shadow_20,type_ZmFuZ3poZW5naGVpdGk=/resize,m_fixed,w_1184)
如上图,如何利用点B及其左右一个点的坐标(共三个点的坐标),来估计点B处一阶导数和二阶导数呢?想法非常简单:既然三个点的坐标已知,那么三个点可以插值得到一个一元二次多项式,求得一元二次多项式的系数后,就可以求多项式在B点处的一阶和二阶导数,从而得到一阶导数和二阶导数的二阶中心差分格式。利用matlab的符号计算可以很容易推导得到。 一阶导数的二阶中心差分格式: ![python二阶中心差法 二阶中心差分推导_二阶导数的四阶中心差分格式_08](https://math-api.51cto.com/?from=%20%20%20%20%20%20%20%20%20f)
二阶导数的二阶中心差分格式: ![python二阶中心差法 二阶中心差分推导_二阶导数的二阶中心差分格式_09](https://math-api.51cto.com/?from=%20%20%20%20%20%20%20%20%20f) 三、一阶导数和二阶导数的四阶中心差分格式![在这里插入图片描述 python二阶中心差法 二阶中心差分推导_一阶导数的二阶中心差分格式_10](https://s2.51cto.com/images/blog/202310/26162531_653a227bdd7d598899.png?x-oss-process=image/watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_30,g_se,x_10,y_10,shadow_20,type_ZmFuZ3poZW5naGVpdGk=/resize,m_fixed,w_1184)
与一阶导数和二阶导数的二阶中心差分推导过程一样,一阶导数和二阶导数的四阶中心差分推导,只需将一元二次多项式插值改为一元四次多项式插值即可。 一阶导数的四阶中心差分格式: ![python二阶中心差法 二阶中心差分推导_二阶导数的四阶中心差分格式_11](https://math-api.51cto.com/?from=%20%20%20%20%20%20%20%20%20f)
二阶导数的四阶中心差分格式: ![python二阶中心差法 二阶中心差分推导_一阶导数的二阶中心差分格式_12](https://math-api.51cto.com/?from=%20%20%20%20%20%20%20%20%20f) 四、MATLAB代码clc
clear
close all
syms x x0 y0 y1 y2 y3 y4 h real
%% 一阶导数和二阶导数的二阶中心差分格式
a = [1, x0, x0^2; 1, (x0 + h), (x0 + h)^2; 1, (x0 + 2 * h), (x0 + 2 * h)^2] \ [y0; y1; y2]; %一元二次多项式y(x) = a1 + a2 * x + a3 * x^2的系数
y(x) = a(1) + a(2) * x + a(3) * x^2;
dy(x) = diff(y, 1);
ddy(x) = diff(y, 2);
dy_two_order_central_difference = simplify(dy(x0 + h))
ddy_two_order_central_difference = simplify(ddy(x0 + h))
%% 一阶导数和二阶导数的四阶中心差分格式
a = [1, x0, x0^2, x0^3, x0^4; 1, (x0 + h), (x0 + h)^2, (x0 + h)^3, (x0 + h)^4; 1, (x0 + 2 * h), (x0 + 2 * h)^2, (x0 + 2 * h)^3, (x0 + 2 * h)^4; ...
1, (x0 + 3 * h), (x0 + 3 * h)^2, (x0 + 3 * h)^3, (x0 + 3 * h)^4; 1, (x0 + 4 * h), (x0 + 4 * h)^2, (x0 + 4 * h)^3, (x0 + 4 * h)^4] \ [y0; y1; y2; y3; y4]; %一元四次多项式y(x) = a1 + a2 * x + a3 * x^2 + a4 * x^3 + a5 * x^4的系数
y(x) = a(1) + a(2) * x + a(3) * x^2 + a(4) * x^3 + a(5) * x^4;
dy(x) = diff(y, 1);
ddy(x) = diff(y, 2);
dy_four_order_central_difference = simplify(dy(x0 + 2 * h))
ddy_four_order_central_difference = simplify(ddy(x0 + 2 * h))
%% 验证
n = 50;
x = linspace(0, 2*pi, n);
h = x(2) - x(1);
y = sin(x);
dy = cos(x);
ddy = -sin(x);
%一阶导数和二阶导数的二阶中心差分(第一个点和最后一个点无法估算一阶导数和二阶导数)
dy1 = nan * zeros(size(x));
ddy1 = nan * zeros(size(x));
for i = 2 : n - 1
dy1(i) = (y(i + 1) - y(i - 1)) / (2.0 * h);
ddy1(i) = (y(i - 1) - 2.0 * y(i) + y(i + 1)) / h^2;
end
%一阶导数和二阶导数的四阶中心差分(第一二个点和最后两个点无法估算一阶导数和二阶导数)
dy2 = nan * zeros(size(x));
ddy2 = nan * zeros(size(x));
for i = 3 : n - 2
dy2(i) = (y(i - 2) - 8.0 * y(i - 1) + 8.0 * y(i + 1) - y(i + 2)) / (12.0 * h);
ddy2(i) = -(y(i - 2) - 16.0 * y(i - 1) + 30.0 * y(i) - 16.0 * y(i + 1) + y(i + 2)) / (12.0 * h^2);
end
max_dy1_err = max(abs(dy1(2 : n - 1) - dy(2 : n - 1)));
max_ddy1_err = max(abs(ddy1(2 : n - 1) - ddy(2 : n - 1)));
max_dy2_err = max(abs(dy2(3 : n - 2) - dy(3 : n - 2)));
max_ddy2_err = max(abs(ddy2(3 : n - 2) - ddy(3 : n - 2)));
disp(['一阶导数的二阶和四阶中心差分近似,最大误差分别为:', num2str(max_dy1_err), ',' , num2str(max_dy2_err)])
disp(['二阶导数的二阶和四阶中心差分近似,最大误差分别为:', num2str(max_ddy1_err), ',' , num2str(max_ddy2_err)])![在这里插入图片描述 python二阶中心差法 二阶中心差分推导_python二阶中心差法_13](https://s2.51cto.com/images/blog/202310/26162532_653a227c08c0270510.png?x-oss-process=image/watermark,size_16,text_QDUxQ1RP5Y2a5a6i,color_FFFFFF,t_30,g_se,x_10,y_10,shadow_20,type_ZmFuZ3poZW5naGVpdGk=/resize,m_fixed,w_1184)
|