【数值分析】二元函数二阶混合偏导数的近似计算式与误差阶推导

二元函数二阶混合偏导数的近似计算式与误差阶推导

问题
引理一：
引理二：
引理三：
引理四：
命题
数值实验
- 函数一
- 函数二
结论

问题

假设f(x,y)f(x,y)f(x,y)在全平面内存在且足够的光滑，求fxy(x0,y0)f_{xy}(x_0,y_0)fxy(x0,y0)的计算式与误差阶

引理一：

若fxy(x0,y0),fyx(x0,y0)f_{xy}(x_0,y_0),f_{yx}(x_0,y_0)fxy(x0,y0),fyx(x0,y0)均在(x0,y0)(x_0,y_0)(x0,y0)处存在且连续，则fxy(x0,y0)=fyx(x0,y0)f_{xy}(x_0,y_0)=f_{yx}(x_0,y_0)fxy(x0,y0)=fyx(x0,y0)

引理二：

设D⊂R2D \subset R^2D⊂R2为一区域，函数f(x,y)∈C2(D)f(x,y)\in C^{2}(D)f(x,y)∈C2(D),且(x0,y0)∈D,(x0+h,y0+k)∈D(x_0,y_0)\in D,(x_0+h,y_0+k) \in D(x0,y0)∈D,(x0+h,y0+k)∈D，则有二元函数带皮亚诺余项的泰勒公式
f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)2+o(ρ2)f(x_0+h,y_0+k)=f(x_0,y_0)+hf_x(x_0,y_0)+kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)+2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2}+o(\rho^2)f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+2h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)+o(ρ2)
其中，ρ=h2+k2\rho=\sqrt{h^2+k^2}ρ=h2+k2

引理三：

(一元函数带拉格朗日余项的泰勒公式)设函数f(x)在包含点x0x_0x0的区间(a,b)内具有直到n+1阶的导数，则对任意x∈(a,b)x\in(a,b)x∈(a,b)都有
f(x)=f(x0)+f′(x0)(x−x0)+12!f′′(x0)(x−x0)2+...+1n!f(n)(x0)(x−x0)n+Rn(x)f(x)=f(x_0)+f^{'}(x_0)(x-x_0)+\frac{1}{2!}f^{''}(x_0)(x-x_0)^2+...+\frac{1}{n!}f^{(n)}(x_0)(x-x_0)^n+R_n(x)f(x)=f(x0)+f′(x0)(x−x0)+2!1f′′(x0)(x−x0)2+...+n!1f(n)(x0)(x−x0)n+Rn(x)
其中,Rn(x)=1(n+1)!f(n+1)(ξ)(x−x0)(n+1)R_n(x)=\frac{1}{(n+1)!}f^{(n+1)}(\xi)(x-x_0)^{(n+1)}Rn(x)=(n+1)!1f(n+1)(ξ)(x−x0)(n+1),ξ\xiξ在x0x_0x0与xxx之间

引理四：

(二元函数带拉格朗日余项的泰勒公式)设D⊂R3D \subset R^3D⊂R3为一区域，函数f(x,y)∈C3(D)f(x,y)\in C^{3}(D)f(x,y)∈C3(D),且(x0,y0)∈D,(x0+h,y0+k)∈D(x_0,y_0)\in D,(x_0+h,y_0+k) \in D(x0,y0)∈D,(x0+h,y0+k)∈D，则有
f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)2!f(x_0+h,y_0+k)=f(x_0,y_0)+hf_x(x_0,y_0)+kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)+2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2!}f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+2!h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)

+h3fxxx(x0+θh,y0+θk)+3h2kfxxy(x0+θh,y0+θk)+3hk2fxyy(x0+θh,y0+θk)+k3fyyy(x0+θh,y0+θk)3!,θ∈(0,1)+\frac{h^3f_{xxx}(x_0+\theta h,y_0+\theta k)+3h^2kf_{xxy}(x_0+\theta h,y_0+\theta k)+3hk^2f_{xyy}(x_0+\theta h,y_0+\theta k)+k^3f_{yyy}(x_0+\theta h,y_0+\theta k)}{3!},\theta \in (0,1)+3!h3fxxx(x0+θh,y0+θk)+3h2kfxxy(x0+θh,y0+θk)+3hk2fxyy(x0+θh,y0+θk)+k3fyyy(x0+θh,y0+θk),θ∈(0,1)

命题

fxy(x0,y0)=f(x0+h,y0+k)−f(x0−h,y0+k)−f(x0+h,y0−k)+f(x0−h,y0−k)4hkf_{xy}(x_0,y_0)=\frac{f(x_0+h,y_0+k)-f(x_0-h,y_0+k)-f(x_0+h,y_0-k)+f(x_0-h,y_0-k)}{4hk}fxy(x0,y0)=4hkf(x0+h,y0+k)−f(x0−h,y0+k)−f(x0+h,y0−k)+f(x0−h,y0−k)

证明：
由引理四得到

f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)2!f(x_0+h,y_0+k)=f(x_0,y_0)+hf_x(x_0,y_0)+kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)+2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2!}f(x0+h,y0+k)=f(x0,y0)+hfx(x0,y0)+kfy(x0,y0)+2!h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)

+h3fxxx(x0+θ1h,y0+θ1k)+3h2kfxxy(x0+θ1h,y0+θ1k)+3hk2fxyy(x0+θ1h,y0+θ1k)+k3fyyy(x0+θ1h,y0+θ1k)3!+\frac{h^3f_{xxx}(x_0+\theta_1 h,y_0+\theta_1 k)+3h^2kf_{xxy}(x_0+\theta_1 h,y_0+\theta_1 k)+3hk^2f_{xyy}(x_0+\theta_1 h,y_0+\theta_1 k)+k^3f_{yyy}(x_0+\theta_1 h,y_0+\theta_1 k)}{3!}+3!h3fxxx(x0+θ1h,y0+θ1k)+3h2kfxxy(x0+θ1h,y0+θ1k)+3hk2fxyy(x0+θ1h,y0+θ1k)+k3fyyy(x0+θ1h,y0+θ1k)

f(x0−h,y0+k)=f(x0,y0)−hfx(x0,y0)+kfy(x0,y0)+h2fxx(x0,y0)−2hkfxy(x0,y0)+k2fyy(x0,y0)2!f(x_0-h,y_0+k)=f(x_0,y_0)-hf_x(x_0,y_0)+kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)-2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2!}f(x0−h,y0+k)=f(x0,y0)−hfx(x0,y0)+kfy(x0,y0)+2!h2fxx(x0,y0)−2hkfxy(x0,y0)+k2fyy(x0,y0)

+−h3fxxx(x0−θ2h,y0+θ2k)+3h2kfxxy(x0−θ2h,y0+θ2k)−3hk2fxyy(x0−θ2h,y0+θ2k)+k3fyyy(x0−θ2h,y0+θ2k)3!+\frac{-h^3f_{xxx}(x_0-\theta_2 h,y_0+\theta_2 k)+3h^2kf_{xxy}(x_0-\theta_2 h,y_0+\theta_2 k)-3hk^2f_{xyy}(x_0-\theta_2 h,y_0+\theta_2 k)+k^3f_{yyy}(x_0-\theta_2 h,y_0+\theta_2 k)}{3!}+3!−h3fxxx(x0−θ2h,y0+θ2k)+3h2kfxxy(x0−θ2h,y0+θ2k)−3hk2fxyy(x0−θ2h,y0+θ2k)+k3fyyy(x0−θ2h,y0+θ2k)

f(x0+h,y0−k)=f(x0,y0)+hfx(x0,y0)−kfy(x0,y0)+h2fxx(x0,y0)−2hkfxy(x0,y0)+k2fyy(x0,y0)2!f(x_0+h,y_0-k)=f(x_0,y_0)+hf_x(x_0,y_0)-kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)-2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2!}f(x0+h,y0−k)=f(x0,y0)+hfx(x0,y0)−kfy(x0,y0)+2!h2fxx(x0,y0)−2hkfxy(x0,y0)+k2fyy(x0,y0)

+h3fxxx(x0+θ3h,y0−θ3k)−3h2kfxxy(x0+θ3h,y0−θ3k)+3hk2fxyy(x0+θ3h,y0−θ3k)−k3fyyy(x0+θ3h,y0−θ3k)3!+\frac{h^3f_{xxx}(x_0+\theta_3 h,y_0-\theta_3 k)-3h^2kf_{xxy}(x_0+\theta_3 h,y_0-\theta_3 k)+3hk^2f_{xyy}(x_0+\theta_3 h,y_0-\theta_3 k)-k^3f_{yyy}(x_0+\theta_3 h,y_0-\theta_3 k)}{3!}+3!h3fxxx(x0+θ3h,y0−θ3k)−3h2kfxxy(x0+θ3h,y0−θ3k)+3hk2fxyy(x0+θ3h,y0−θ3k)−k3fyyy(x0+θ3h,y0−θ3k)

f(x0−h,y0−k)=f(x0,y0)−hfx(x0,y0)−kfy(x0,y0)+h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)2!f(x_0-h,y_0-k)=f(x_0,y_0)-hf_x(x_0,y_0)-kf_y(x_0,y_0)+\frac{h^2f_{xx}(x_0,y_0)+2hkf_{xy}(x_0,y_0)+k^2f_{yy}(x_0,y_0)}{2!}f(x0−h,y0−k)=f(x0,y0)−hfx(x0,y0)−kfy(x0,y0)+2!h2fxx(x0,y0)+2hkfxy(x0,y0)+k2fyy(x0,y0)

+−h3fxxx(x0−θ4h,y0−θ4k)−3h2kfxxy(x0−θ4h,y0−θ4k)−3hk2fxyy(x0−θ4h,y0−θ4k)−k3fyyy(x0−θ4h,y0−θ4k)3!+\frac{-h^3f_{xxx}(x_0-\theta_4 h,y_0-\theta_4 k)-3h^2kf_{xxy}(x_0-\theta_4 h,y_0-\theta_4 k)-3hk^2f_{xyy}(x_0-\theta_4 h,y_0-\theta_4 k)-k^3f_{yyy}(x_0-\theta_4 h,y_0-\theta_4 k)}{3!}+3!−h3fxxx(x0−θ4h,y0−θ4k)−3h2kfxxy(x0−θ4h,y0−θ4k)−3hk2fxyy(x0−θ4h,y0−θ4k)−k3fyyy(x0−θ4h,y0−θ4k)

所以，误差等于
f(x0+h,y0+k)−f(x0+h,y0−k)−f(x0−h,y0+k)+f(x0−h,y0−k)4hk−fyx(x0,y0)\frac{f(x_0+h,y_0+k)-f(x_0+h,y_0-k)-f(x_0-h,y_0+k)+f(x_0-h,y_0-k)}{4hk}-f_{yx}(x_0,y_0)4hkf(x0+h,y0+k)−f(x0+h,y0−k)−f(x0−h,y0+k)+f(x0−h,y0−k)−fyx(x0,y0)

=124hk((h3fxxx(x0+θ1h,y0+θ1k)+3h2kfxxy(x0+θ1h,y0+θ1k)+3hk2fxyy(x0+θ1h,y0+θ1k)+k3fyyy(x0+θ1h,y0+θ1k))=\frac{1}{24hk}((h^3f_{xxx}(x_0+\theta_1 h,y_0+\theta_1 k)+3h^2kf_{xxy}(x_0+\theta_1 h,y_0+\theta_1 k)+3hk^2f_{xyy}(x_0+\theta_1 h,y_0+\theta_1 k)+k^3f_{yyy}(x_0+\theta_1 h,y_0+\theta_1 k))=24hk1((h3fxxx(x0+θ1h,y0+θ1k)+3h2kfxxy(x0+θ1h,y0+θ1k)+3hk2fxyy(x0+θ1h,y0+θ1k)+k3fyyy(x0+θ1h,y0+θ1k))

−(−h3fxxx(x0−θ2h,y0+θ2k)+3h2kfxxy(x0−θ2h,y0+θ2k)−3hk2fxyy(x0−θ2h,y0+θ2k)+k3fyyy(x0−θ2h,y0+θ2k))-(-h^3f_{xxx}(x_0-\theta_2 h,y_0+\theta_2 k)+3h^2kf_{xxy}(x_0-\theta_2 h,y_0+\theta_2 k)-3hk^2f_{xyy}(x_0-\theta_2 h,y_0+\theta_2 k)+k^3f_{yyy}(x_0-\theta_2 h,y_0+\theta_2 k))−(−h3fxxx(x0−θ2h,y0+θ2k)+3h2kfxxy(x0−θ2h,y0+θ2k)−3hk2fxyy(x0−θ2h,y0+θ2k)+k3fyyy(x0−θ2h,y0+θ2k))

−(h3fxxx(x0+θ3h,y0−θ3k)−3h2kfxxy(x0+θ3h,y0−θ3k)+3hk2fxyy(x0+θ3h,y0−θ3k)−k3fyyy(x0+θ3h,y0−θ3k))-(h^3f_{xxx}(x_0+\theta_3 h,y_0-\theta_3 k)-3h^2kf_{xxy}(x_0+\theta_3 h,y_0-\theta_3 k)+3hk^2f_{xyy}(x_0+\theta_3 h,y_0-\theta_3 k)-k^3f_{yyy}(x_0+\theta_3 h,y_0-\theta_3 k))−(h3fxxx(x0+θ3h,y0−θ3k)−3h2kfxxy(x0+θ3h,y0−θ3k)+3hk2fxyy(x0+θ3h,y0−θ3k)−k3fyyy(x0+θ3h,y0−θ3k))

+(−h3fxxx(x0−θ4h,y0−θ4k)−3h2kfxxy(x0−θ4h,y0−θ4k)−3hk2fxyy(x0−θ4h,y0−θ4k)−k3fyyy(x0−θ4h,y0−θ4k)))+(-h^3f_{xxx}(x_0-\theta_4 h,y_0-\theta_4 k)-3h^2kf_{xxy}(x_0-\theta_4 h,y_0-\theta_4 k)-3hk^2f_{xyy}(x_0-\theta_4 h,y_0-\theta_4 k)-k^3f_{yyy}(x_0-\theta_4 h,y_0-\theta_4 k)))+(−h3fxxx(x0−θ4h,y0−θ4k)−3h2kfxxy(x0−θ4h,y0−θ4k)−3hk2fxyy(x0−θ4h,y0−θ4k)−k3fyyy(x0−θ4h,y0−θ4k)))

将等式右边合并同类项按照h3,h2k,hk2,k3h^3,h^2k,hk^2,k^3h3,h2k,hk2,k3分成四类考虑

含h3h^3h3的几项为

h324hk(fxxx(x0+θ1h,y0+θ1k)+fxxx(x0−θ2h,y0+θ2k)−fxxx(x0+θ3h,y0−θ3k)−fxxx(x0−θ4h,y0−θ4k))\frac{h^3}{24hk}(f_{xxx}(x_0+\theta_1 h,y_0+\theta_1 k)+f_{xxx}(x_0-\theta_2 h,y_0+\theta_2 k)-f_{xxx}(x_0+\theta_3 h,y_0-\theta_3 k)-f_{xxx}(x_0-\theta_4 h,y_0-\theta_4 k))24hkh3(fxxx(x0+θ1h,y0+θ1k)+fxxx(x0−θ2h,y0+θ2k)−fxxx(x0+θ3h,y0−θ3k)−fxxx(x0−θ4h,y0−θ4k))

设G(θ)=fxxx(x0+θh,y0+θk),H(θ)=fxxx(x0−θh,y0+θk)G(\theta)=f_{xxx}(x_0+\theta h,y_0+\theta k),H(\theta)=f_{xxx}(x_0-\theta h,y_0+\theta k)G(θ)=fxxx(x0+θh,y0+θk),H(θ)=fxxx(x0−θh,y0+θk)

则上面几项可以化为h224k((G(θ1)−G(−θ4)+(H(θ2)−H(−θ3)))\frac{h^2}{24k}((G(\theta_1)-G(-\theta_4)+(H(\theta_2)-H(-\theta_3)))24kh2((G(θ1)−G(−θ4)+(H(θ2)−H(−θ3)))

=h224k((θ1+θ4)G′(ξ1)+(θ2+θ3)H′(ξ2))=\frac{h^2}{24k}((\theta_1+\theta_4)G^{'}(\xi_1)+(\theta_2+\theta_3)H^{'}(\xi_2))=24kh2((θ1+θ4)G′(ξ1)+(θ2+θ3)H′(ξ2))

而用链式法则对G,HG,HG,H求导，得到

G′(θ)=hfxxxx(x0+θh,y0+θk)+kfxxxy(x0+θh,y0+θk)G^{'}(\theta)=hf_{xxxx}(x_0+\theta h,y_0+\theta k)+kf_{xxxy}(x_0+\theta h,y_0+\theta k)G′(θ)=hfxxxx(x0+θh,y0+θk)+kfxxxy(x0+θh,y0+θk)

H′(θ)=−hfxxxx(x0−θh,y0+θk)+kfxxxy(x0−θh,y0+θk)H^{'}(\theta)=-hf_{xxxx}(x_0-\theta h,y_0+\theta k)+kf_{xxxy}(x_0-\theta h,y_0+\theta k)H′(θ)=−hfxxxx(x0−θh,y0+θk)+kfxxxy(x0−θh,y0+θk)

因此，原来几项可以化为

h224k((θ1+θ4)(hfxxxx(x0+ξ1h,y0+ξ1k)+kfxxxy(x0+ξ1h,y0+ξ1k))\frac{h^2}{24k}((\theta_1+\theta_4)(hf_{xxxx}(x_0+\xi_1 h,y_0+\xi_1 k)+kf_{xxxy}(x_0+\xi_1 h,y_0+\xi_1 k))24kh2((θ1+θ4)(hfxxxx(x0+ξ1h,y0+ξ1k)+kfxxxy(x0+ξ1h,y0+ξ1k))

+(θ2+θ3)(−hfxxxx(x0−ξ2h,y0+ξ2k)+kfxxxy(x0−ξ2h,y0+ξ2k)))+(\theta_2+\theta_3)(-hf_{xxxx}(x_0-\xi_2 h,y_0+\xi_2 k)+kf_{xxxy}(x_0-\xi_2 h,y_0+\xi_2 k)))+(θ2+θ3)(−hfxxxx(x0−ξ2h,y0+ξ2k)+kfxxxy(x0−ξ2h,y0+ξ2k)))

由于θ1+θ4\theta_1+\theta_4θ1+θ4与θ2+θ3\theta_2+\theta_3θ2+θ3一般不相等，所以无法用拉格朗日中值定理进一步化简

因此当h,kh,kh,k大小接近的时候，若假设上面的四阶偏导数均有界，则有误差阶数为O(h2)O(h^2)O(h2)

类似地对h2k,hk2,k3h^2k,hk^2,k^3h2k,hk2,k3讨论，得到在h与k同阶的时候，误差阶数为O(h2)O(h^2)O(h2)

数值实验

函数一

f(x,y)=x2sin(y)+exy2f(x,y)=x^2sin(y)+e^{xy^2}f(x,y)=x2sin(y)+exy2

fx(x,y)=2xsin(y)+y2exy2f_{x}(x,y)=2xsin(y)+y^2e^{xy^2}fx(x,y)=2xsin(y)+y2exy2

fy(x,y)=x2cos(y)+2xyexy2f_{y}(x,y)=x^2cos(y)+2xye^{xy^2}fy(x,y)=x2cos(y)+2xyexy2

fxy(x,y)=2xcos(y)+2yexy2+2xy3exy2f_{xy}(x,y)=2xcos(y)+2ye^{xy^2}+2xy^3e^{xy^2}fxy(x,y)=2xcos(y)+2yexy2+2xy3exy2

 import numpy as npdef f1(x,y):return x*x*np.sin(y)+np.exp(x*y*y)def f1xy(x,y):return 2*x*np.cos(y)+2*y*np.exp(x*y*y)+2*x*y*y*y*np.exp(x*y*y)x=2.5y=1.8h=0.1k=0.1print("(%.4f,%.4f)二阶混合偏导数真实值:%.10f"%(x,y,f2xy(x,y)))print("序号\th,k\t\t近似公式计算值\t\t误差\t\t误差衰减倍数")for i in range(20):now=(f1(x+h,y+k)-f1(x-h,y+k)-f1(x+h,y-k)+f1(x-h,y-k))*0.25/h/kprint("%d\t%.7f\t%.10f\t%.10f"%(i+1,h,now,now-f1xy(x,y)),end='')if i>0:print("\t%.6f"%(las/(now-f1xy(x,y))))else:print()las=now-f1xy(x,y)h*=0.5k*=0.5

(2.5000,1.8000)二阶混合偏导数真实值:2.6708851424
序号 h,k 近似公式计算值误差误差衰减倍数
1 0.1000000 133118.5135877089 25192.8754518838
2 0.0500000 113820.5496154785 5894.9114796533 4.273665
3 0.0250000 109375.4104465985 1449.7723107733 4.066095
4 0.0125000 108286.6028087446 360.9646729195 4.016383
5 0.0062500 108015.7871949166 90.1490590915 4.004087
6 0.0031250 107948.1696480652 22.5315122401 4.001021
7 0.0015625 107931.2706538942 5.6325180690 4.000256
8 0.0007813 107927.0462427288 1.4081069037 4.000064
9 0.0003906 107925.9901575744 0.3520217493 4.000057
10 0.0001953 107925.7261663675 0.0880305424 3.998859
11 0.0000977 107925.6601095200 0.0219736948 4.006178
12 0.0000488 107925.6436347961 0.0054989710 3.995965
13 0.0000244 107925.6391525269 0.0010167017 5.408637
14 0.0000122 107925.6431579590 0.0050221338 0.202444
15 0.0000061 107925.5615234375 -0.0766123876 -0.065553
16 0.0000031 107925.3906250000 -0.2475108251 0.309531
17 0.0000015 107925.4882812500 -0.1498545751 1.651673
18 0.0000008 107926.7578125000 1.1196766749 -0.133837
19 0.0000004 107917.9687500000 -7.6693858251 -0.145993
20 0.0000002 107887.5000000000 -38.1381358251 0.201095

函数二

f(x,y)=x32sin(x+y)f(x,y)=x^{\frac{3}{2}}sin(x+y)f(x,y)=x23sin(x+y)

fx(x,y)=x32cos(x+y)+32x12sin(x+y)f_{x}(x,y)=x^{\frac{3}{2}}cos(x+y)+\frac{3}{2}x^{\frac{1}{2}}sin(x+y)fx(x,y)=x23cos(x+y)+23x21sin(x+y)

fy(x,y)=x32cos(x+y)f_{y}(x,y)=x^{\frac{3}{2}}cos(x+y)fy(x,y)=x23cos(x+y)

fxy(x,y)=−x32sin(x+y)+32x12cos(x+y)f_{xy}(x,y)=-x^{\frac{3}{2}}sin(x+y)+\frac{3}{2}x^{\frac{1}{2}}cos(x+y)fxy(x,y)=−x23sin(x+y)+23x21cos(x+y)

 import numpy as npdef f2(x,y):return x**1.5*np.sin(x+y)def f2xy(x,y):return 1.5*(x**0.5)*np.cos(x+y)-x**1.5*np.sin(x+y)x=5y=5h=0.1k=0.1print("(%.4f,%.4f)二阶混合偏导数真实值:%.10f"%(x,y,f2xy(x,y)))print("序号\th,k\t\t近似公式计算值\t\t误差\t\t误差衰减倍数")for i in range(15):now=(f2(x+h,y+k)-f2(x-h,y+k)-f2(x+h,y-k)+f2(x-h,y-k))*0.25/h/kprint("%d\t%.7f\t%.10f\t%.10f"%(i+1,h,now,now-f2xy(x,y)),end='')if i>0:print("\t%.6f"%(las/(now-f2xy(x,y))))else:print()las=now-f2xy(x,y)h*=0.5k*=0.5

(5.0000,5.0000)二阶混合偏导数真实值:3.2680094602
序号 h,k 近似公式计算值误差误差衰减倍数
1 0.1000000 3.2674426586 -0.0005668017
2 0.0500000 3.2678703460 -0.0001391142 4.074362
3 0.0250000 3.2679748435 -0.0000346167 4.018704
4 0.0125000 3.2680008162 -0.0000086441 4.004682
5 0.0062500 3.2680072998 -0.0000021604 4.001154
6 0.0031250 3.2680089201 -0.0000005401 4.000008
7 0.0015625 3.2680093253 -0.0000001350 4.001843
8 0.0007813 3.2680094264 -0.0000000338 3.989878
9 0.0003906 3.2680094504 -0.0000000098 3.446221
10 0.0001953 3.2680094533 -0.0000000069 1.421488
11 0.0000977 3.2680094242 -0.0000000360 0.191759
12 0.0000488 3.2680093311 -0.0000001291 0.278833
13 0.0000244 3.2680094242 -0.0000000360 3.586371
14 0.0000122 3.2680064440 -0.0000030162 0.011938
15 0.0000061 3.2679975033 -0.0000119569 0.252259

由上面几个函数的例子可以看出，h和k减小为原来一半时，误差减小为原来的四分之一，因此验证了误差的阶数为平方级别

结论

对于足够光滑的函数f(x,y)f(x,y)f(x,y),有二阶混合偏导数计算式：

该计算式的误差阶上界数为二阶

由于作者水平有限，如果推导过程中有错误或考虑不周之处，还望不吝指正。