问题

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

引理一：

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

引理二：

设$D\subset R^2$为一区域，函数$f(x,y)\in C^2(D)$,且$(x_0,y_0)\in D,(x_0+h,y_0+k)\in D$，则有二元函数带皮亚诺余项的泰勒公式
$$f(x_0+h,y_0+k)=f(x_0,y_0)+h{f}_x(x_0,y_0)+k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)+2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2}+o({\rho }^2)$$
其中，$\rho =\sqrt{h^2+k^2}$

引理三：

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

引理四：

(二元函数带拉格朗日余项的泰勒公式)设$D\subset R^3$为一区域，函数$f(x,y)\in C^3(D)$,且$(x_0,y_0)\in D,(x_0+h,y_0+k)\in D$，则有
$f(x_0+h,y_0+k)=f(x_0,y_0)+h{f}_x(x_0,y_0)+k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)+2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2!}$

$+\frac{h^3{f}_{xxx}(x_0+\theta h,y_0+\theta k)+3h^{2}k{f}_{xxy}(x_0+\theta h,y_0+\theta k)+3h{k}^2{f}_{xyy}(x_0+\theta h,y_0+\theta k)+k^3{f}_{yyy}(x_0+\theta h,y_0+\theta k)}{3!},\theta \in (0,1)$

命题

$f_{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}$

证明：
由引理四得到

$f(x_0+h,y_0+k)=f(x_0,y_0)+h{f}_x(x_0,y_0)+k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)+2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2!}$

$+\frac{h^3{f}_{xxx}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+3h^{2}k{f}_{xxy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+3h{k}^2{f}_{xyy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+k^3{f}_{yyy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)}{3!}$

$f(x_0-h,y_0+k)=f(x_0,y_0)-h{f}_x(x_0,y_0)+k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)-2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2!}$

$+\frac{-h^3{f}_{xxx}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)+3h^{2}k{f}_{xxy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)-3h{k}^2{f}_{xyy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)+k^3{f}_{yyy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)}{3!}$

$f(x_0+h,y_0-k)=f(x_0,y_0)+h{f}_x(x_0,y_0)-k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)-2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2!}$

$+\frac{h^3{f}_{xxx}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)-3h^{2}k{f}_{xxy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)+3h{k}^2{f}_{xyy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)-k^3{f}_{yyy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)}{3!}$

$f(x_0-h,y_0-k)=f(x_0,y_0)-h{f}_x(x_0,y_0)-k{f}_y(x_0,y_0)+\frac{h^2{f}_{xx}(x_0,y_0)+2hk{f}_{xy}(x_0,y_0)+k^2{f}_{yy}(x_0,y_0)}{2!}$

$+\frac{-h^3{f}_{xxx}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-3h^{2}k{f}_{xxy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-3h{k}^2{f}_{xyy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-k^3{f}_{yyy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)}{3!}$

所以，误差等于
$\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)$

$=\frac{1}{24hk}((h^3{f}_{xxx}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+3h^{2}k{f}_{xxy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+3h{k}^2{f}_{xyy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k)+k^3{f}_{yyy}(x_0+{\theta }_{1}h,y_0+{\theta }_{1}k))$

$-(-h^3{f}_{xxx}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)+3h^{2}k{f}_{xxy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)-3h{k}^2{f}_{xyy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k)+k^3{f}_{yyy}(x_0-{\theta }_{2}h,y_0+{\theta }_{2}k))$

$-(h^3{f}_{xxx}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)-3h^{2}k{f}_{xxy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)+3h{k}^2{f}_{xyy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k)-k^3{f}_{yyy}(x_0+{\theta }_{3}h,y_0-{\theta }_{3}k))$

$+(-h^3{f}_{xxx}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-3h^{2}k{f}_{xxy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-3h{k}^2{f}_{xyy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)-k^3{f}_{yyy}(x_0-{\theta }_{4}h,y_0-{\theta }_{4}k)))$

将等式右边合并同类项按照$h^3,h^{2}k,h{k}^2,k^3$分成四类考虑

含$h^3$的几项为

$\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))$

设$G(\theta )=f_{xxx}(x_0+\theta h,y_0+\theta k),H(\theta )=f_{xxx}(x_0-\theta h,y_0+\theta k)$

则上面几项可以化为$\frac{h^2}{24k}((G({\theta }_1)-G(-{\theta }_4)+(H({\theta }_2)-H(-{\theta }_3)))$

$=\frac{h^2}{24k}(({\theta }_1+{\theta }_4)G^{{}^{\prime }}({\xi }_1)+({\theta }_2+{\theta }_3)H^{{}^{\prime }}({\xi }_2))$

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

$G^{{}^{\prime }}(\theta )=h{f}_{xxxx}(x_0+\theta h,y_0+\theta k)+k{f}_{xxxy}(x_0+\theta h,y_0+\theta k)$

$H^{{}^{\prime }}(\theta )=-h{f}_{xxxx}(x_0-\theta h,y_0+\theta k)+k{f}_{xxxy}(x_0-\theta h,y_0+\theta k)$

因此，原来几项可以化为

$\frac{h^2}{24k}(({\theta }_1+{\theta }_4)(h{f}_{xxxx}(x_0+{\xi }_{1}h,y_0+{\xi }_{1}k)+k{f}_{xxxy}(x_0+{\xi }_{1}h,y_0+{\xi }_{1}k))$

$+({\theta }_2+{\theta }_3)(-h{f}_{xxxx}(x_0-{\xi }_{2}h,y_0+{\xi }_{2}k)+k{f}_{xxxy}(x_0-{\xi }_{2}h,y_0+{\xi }_{2}k)))$

由于${\theta }_1+{\theta }_4$与${\theta }_2+{\theta }_3$一般不相等，所以无法用拉格朗日中值定理进一步化简

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

类似地对$h^{2}k,h{k}^2,k^3$讨论，得到在h与k同阶的时候，误差阶数为$O(h^2)$

数值实验

函数一

$f(x,y)=x^{2}sin(y)+e^{x{y}^2}$

$f_x(x,y)=2xsin(y)+y^2{e}^{x{y}^2}$

$f_y(x,y)=x^{2}cos(y)+2xy{e}^{x{y}^2}$

$f_{xy}(x,y)=2xcos(y)+2y{e}^{x{y}^2}+2x{y}^3{e}^{x{y}^2}$

	import numpy as np
	def 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.5
	y=1.8
	h=0.1
	k=0.1
	print("(%.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/k
	    print("%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.5
	    k*=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)=x^{\frac{3}{2}}sin(x+y)$

$f_x(x,y)=x^{\frac{3}{2}}cos(x+y)+\frac{3}{2}{x}^{\frac{1}{2}}sin(x+y)$

$f_y(x,y)=x^{\frac{3}{2}}cos(x+y)$

$f_{xy}(x,y)=-x^{\frac{3}{2}}sin(x+y)+\frac{3}{2}{x}^{\frac{1}{2}}cos(x+y)$

	import numpy as np
	def 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=5
	y=5
	h=0.1
	k=0.1
	print("(%.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/k
	    print("%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.5
	    k*=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_{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}$$

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

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