kill(all)$

mixed partial derivatives ∂/∂x(∂u/∂y) (numerical approximation)
peter.vlasschaert@gmail.com,03/11/2017

to find ?  mixed partial derivatives ∂/∂x(∂u/∂y) 
suppose : polynomial (x,y)    'y(x,y)'  => (x,y,x*y,x^2,y^2)

y(x,y):=a_[1]*x+a_[2]*y+a_[3]*x*y+a_[4]*x^2+a_[5]*y^2;

rem : we only need 5 equations 
   general :polynomial (x,y)  => (1,x,y,x*y,x^2,y^2)
rem : general 6 equations we only need 5 equations.

suppose : gridspacing 'x direction ' = h
                    gridspacing 'y direction ' = k

p1:y1=y(x0+h,y0+k);
p2:y2=y(x0+h,y0-k);
p3:y3=y(x0-h,y0+k);
p4:y4=y(x0-h,y0-k);
p5:y5=y(x0,y0);

p6:solve([p1,p2,p3,p4,p5], [a_[1],a_[2],a_[3],a_[4],a_[5]])[1]$
p7a:p6[1];
p7b:p6[2];
p7c:p6[3];
p7d:p6[4];

solution from the above: a_[3]
----------------------------------------------
reason : independent of x0,y0   (used for every point  inside the domain ' regular points').

p8a:'diff(u,x,1,y,1) = rhs(p7c);
p8b:'diff(u,y,1,x,1) = rhs(p7c);

p9a:y1=u(x0+h,y0+k);
p9b:y2=u(x0+h,y0-k);
p9c:y3=u(x0-h,y0+k);
p9d:y4=u(x0-h,y0-k);
p9e:y5=u(x0,y0);

numerical representation :∂/∂x(∂u/∂y) 

p10:ev(p8a,p9a,p9b,p9c,p9d,p9e);

simplify when , h=k
result : same grid spacing ' x ' = same grid spacing ' y ' = h

p11:ev(p10,k=h);
p12:ev(p11,x0=i,y0=j);

index representation i,j

uu(i,j):=uu[i,j];

p13:ev(p12,u=uu)$
p14:ev(p13,uu=u);

when h=1 , "  grid spacing with one unit  "

p15:ev(p14,h=1);

rem: finite differences for ∂/∂x(∂u/∂y) 