Partial derivatives & Vector calculus

Partial derivatives

Functions of several arguments (multivariate functions) such as f[x,y] can be differentiated with respect to each argument

((∂f)/(∂x))≡(∂_x)f,        ((∂f)/(∂y))≡(∂_y)f,

etc. One can define higher-order derivatives with respect to the same or different variables

(((∂^2)f)/(∂(x^2)))≡(∂_(x,x))f,        (((∂^2)f)/(∂(y^2)))≡(∂_(y,y))f,     (((∂^2)f)/(∂x∂y))≡(∂/(∂x))((∂f)/(∂y))≡(∂_(x,y))f

For most of the functions mixed partial derivatives do not depend on the order of differentiation

(((∂^2)f)/(∂x∂y))≡(((∂^2)f)/(∂y∂x))

This holds if the mixed derivatives are continuous at a given point. For instance,

f[x_,y_]=x y;
(∂_x)f[x,y]
(∂_y)f[x,y]
(∂_(y,x))f[x,y]
(∂_(x,y))f[x,y]

"Bad" functions

An example of a "bad" function is

f[x_,y_]=x y (((x^2)-(y^2))/((x^2)+(y^2)));
fx[x_,y_]=Simplify[(∂_x)f[x,y]]
fy[x_,y_]=Simplify[(∂_y)f[x,y]]
fxy[x_,y_]=Simplify[(∂_(y,x))f[x,y]]
fyx[x_,y_]=Simplify[(∂_(x,y))f[x,y]]

Here the mixed derivative do not depend on the order everywhere except (0,0). At (0,0) they do not have well defined values and depend on the way (x,y) approaches to (0,0)

Limit[fxy[x,y],x->0]
Limit[fxy[x,y],y->0]

It is convenient to rewrite these functions in the polar coordinates (r,ϕ) to elucidate their behavior at r=0

ff[r_,ϕ_]=Simplify[f[x,y]/.{x->r Cos[ϕ],y->r Sin[ϕ]}]

This has a well-defined limit r->0. Further

ffx[r_,ϕ_]=Simplify[fx[x,y]/.{x->r Cos[ϕ],y->r Sin[ϕ]}]
ffy[r_,ϕ_]=Simplify[fy[x,y]/.{x->r Cos[ϕ],y->r Sin[ϕ]}]

These also tends to zero as  r->0. Finally

ffxy[r_,ϕ_]=Simplify[fxy[x,y]/.{x->r Cos[ϕ],y->r Sin[ϕ]}]

The mixed derivative does not have a well defined value at  r->0 and depends on the ϕ.

Plot3D[f[x,y],{x,-1,1},{y,-1,1},PlotPoints->100]

Plot3D[fx[x,y],{x,-1,1},{y,-1,1},PlotPoints->100]

Plot3D[fxy[x,y],{x,-1,1},{y,-1,1},PlotPoints->100]

Multivariate series

Taylor series can be generalized for multivariate functions and the

f[x_,y_]=Sin[x+y];
Series[f[x,y],{x,0,3}]

Series[f[x,y],{x,0,3},{y,0,3}]

or, in the symmetric form

Normal[Series[f[x,y],{x,0,3},{y,0,3}]]
Expand[Normal[Series[f[x,y],{x,0,3},{y,0,3}]]]

Exercise: Find a way to sort this polynomial in increasing powers of x, y.

Vector calculus

Physics makes use of vector differential operations on functions such as gradient, divergence, curl (rotor), Laplacian, etc. 
In the current version of  Mathematica realizations of these operations are new and not included in the main body of the software. Instead, these functions are implemented in the optional VectorAnalysis package that has to be called before performing these operations

Needs["VectorAnalysis`"]

Unfortunately, this package seems to be inconvenient.

Gradient

Gradient of a scalar funtion is a vector defined as

grad f≡∇f≡(Subscript[e, x])((∂f)/(∂x))+(Subscript[e, y])((∂f)/(∂y))+(Subscript[e, z])((∂f)/(∂z))≡{(∂_x)f,(∂_y)f,(∂_z)f}

One can speak about the gradient operator defined as

∇≡(Subscript[e, x])(∂/(∂x))+(Subscript[e, y])(∂/(∂y))+(Subscript[e, z])(∂/(∂z))

that acts on scalar functions of vector arguments. An example of a gradient in physics is force F that is minus gradient of the potential energy U[x,y,z] and similar for the electric field E that is minus gradient of the electric potential ϕ

F≡-∇U,         E≡-∇ϕ

Examples trying to use the Mathematica's VectorCalculus package:

Following Mathematica help:

Clear[x,y,z,U]
U=(x^2)+(y^2)+(z^2);
Grad[U]

- a wrong output. An attempt of a standard usage

U[x_,y_,z_]=(x^2)+(y^2)+(z^2);   (* 3d oscillator *)

F[x_,y_,z_]=Grad[U[x,y,z]]

- same wrong result. Still this command is working with a special naming choice

Grad[(Xx^2)+(Yy^2)+(Zz^2)]

However, this naming restriction is inconvenient.

Fortunately, it is not difficult to program the gradient in Mathematica. To use the definition below, quit the kernel to remove the VectorAnalysis package from the memory

Grad[f_]:={(∂_x)f,(∂_y)f,(∂_z)f}

Since x,y,z enters the definition of this function, the arguments of f should also be x,y,z. With any other notation for the arguments of f, it won't work, in contrast to definitions of tru functions.
MyGrad works on expressions

Grad[(x^2)+(y^2)+(z^2)]     (* ∇(r^2)=2r *)

This means

∇(r^2)=2r

MyGrad also works on functions

U[x_,y_,z_]=(x^2)+(y^2)+(z^2);     (* 3d oscillator *)
V[x_,y_,z_]=(x^4);   
Grad[U[x,y,z]]
Grad[V[x,y,z]]

Mathematica has the symbol ∇ but it seems it is only for typing

?∇

?E

I==I

One can define a vector function that is the gradient of a scalar function. For the electric field of a point charge one has

ϕ[x_,y_,z_]=((k Q)/(Sqrt[(x^2)+(y^2)+(z^2)]));   (* Coulomb potential of a point charge *)
EE[x_,y_,z_]=-Grad[ϕ[x,y,z]]   (* Electric field of a point charge *)

Divergence

Divergence of a vector is a scalar defined by

div A ≡∇·A≡((∂(Subscript[A, x]))/(∂x))+((∂(Subscript[A, y]))/(∂y))+((∂(Subscript[A, z]))/(∂z))

divergence can be represented by the operator

∇≡(Subscript[e, x])(∂/(∂x))+(Subscript[e, y])(∂/(∂y))+(Subscript[e, z])(∂/(∂z))

same as the gradient operator above. The only difference between them is that gradient acts on scalars and divergence acts on vectors.

Div[Avec_]:=(∂_x)Avec[[1]]+(∂_y)Avec[[2]]+(∂_z)Avec[[3]]

Examples

Div[{x,y,z}]

Fvec[x_,y_,z_]={(x^2),(y^2),(z^2)};
Div[Fvec[x,y,z]]

Fvec[x_,y_,z_]={y,x,x y};
Div[Fvec[x,y,z]]

Laplacian

Laplacian is a second-order vector differential operation. Laplacian of a scalar f is defined as div grad f and denoted by Δ or (∇^2)

Δ f≡div grad f ≡∇·∇ f≡(∇^2)f

From this definition follows

Δ f=(((∂^2)f)/(∂(x^2)))+(((∂^2)f)/(∂(y^2)))+(((∂^2)f)/(∂(z^2)))

The Laplace operator

Δ ≡(∇^2)=((∂^2)/(∂(x^2)))+((∂^2)/(∂(y^2)))+((∂^2)/(∂(z^2)))

can be obtained by squaring the gradient / divergence operator above.

Laplace[f_]:=(∂_(x,x))f+(∂_(y,y))f+(∂_(z,z))f

Example

f[x_,y_,z_]=(x^2)+(y^2)+(z^2);
Laplace[f[x,y,z]]

Laplacian of a vector is defined in components through Laplacians of the components as

ΔA ≡(Subscript[e, 1])(Subscript[ΔA, 1])+ (Subscript[e, 2])(Subscript[ΔA, 2])+(Subscript[e, 3])(Subscript[ΔA, 3])

LaplaceVec[Avec_]:={(∂_(x,x))Avec[[1]]+(∂_(y,y))Avec[[1]]+(∂_(z,z))Avec[[1]],(∂_(x,x))Avec[[2]]+(∂_(y,y))Avec[[2]]+(∂_(z,z))Avec[[2]],(∂_(x,x))Avec[[3]]+(∂_(y,y))Avec[[3]]+(∂_(z,z))Avec[[3]]}

Rotor (curl)

Curl or rotor of a vector is a vector defined by

rot A≡curl A≡det[(({{(Subscript[e, x]), (Subscript[e, y]), (Subscript[e, z])}, {(∂_x), (∂_y), (∂_z)}, {(Subscript[A, x]), (Subscript[A, y]), (Subscript[A, z])}}))]≡∇×A

Rot[Avec_]:={(∂_y)Avec[[3]]-(∂_z)Avec[[2]],(∂_z)Avec[[1]]-(∂_x)Avec[[3]],(∂_x)Avec[[2]]-(∂_y)Avec[[1]]}
Curl[Avec_]:=Rot[Avec]

Examples

Curl[{x,y,z}]

Rot[{x,y,z}]

Curl[{y,-x,0}]

Repeated first-order differential vector operations

Repeated vector differential operations satisfy some identities that are similar to repeated vector products if one uses the  ∇ operator and considers it as a vector.

Double curl

Similarly to the double vector product A×(B×C) ≡ B (A·C) - C (A·B), one has

rot rot A≡∇×(∇×A)≡∇(∇·A)-(∇^2)A≡grad(div A)-Δ A

Note that here the vector A is always rightmost because the differential operators are acting on it. Check this identity

Avec[x_,y_,z_]={Avecx[x,y,z],Avecy[x,y,z],Avecz[x,y,z]};

Rot[Rot[Avec[x,y,z]]]

Grad[Div[Avec[x,y,z]]]

LaplaceVec[Avec[x,y,z]]

Rot[Rot[Avec[x,y,z]]]==Grad[Div[Avec[x,y,z]]]-LaplaceVec[Avec[x,y,z]]

Divergence of a curl

(of a solenoidal field) is zero

div curl A≡∇·(∇×A)≡0

Check this identity

Div[Curl[Avec[x,y,z]]]

Curl of a gradient

(of a potential field) is also zero

rot grad f≡∇×(∇f)≡0

Check this identity

Rot[Grad[f]]

Gradient of a divergence

This one seems to be not expressible via other operations

grad div A≡∇(∇·A)≡?

Grad[Div[Avec[x,y,z]]]

Differential vector operations on two objects

Gradient of a product of two scalars

grad (ϕ ψ)≡∇(ϕ ψ)≡ϕ∇ψ+ψ∇ϕ

Check this identity

Grad[ϕ[x,y,z]ψ[x,y,z]]==ϕ[x,y,z]Grad[ψ[x,y,z]]+Grad[ϕ[x,y,z]]ψ[x,y,z]

Divergence and curl of a product of a vector and a scalar

div (ϕ A)≡∇·(ϕ A)≡ϕ∇·A+∇ϕ·A
rot (ϕ A)≡∇×(ϕ A)≡ϕ∇×A+∇ϕ×A

Excercise: Check these identities

Divergence of a vector product

div (A×B)≡∇·(A×B)≡(∇×A)·B-A·(∇×B)

Check this identity

Avec[x_,y_,z_]={Avecx[x,y,z],Avecy[x,y,z],Avecz[x,y,z]};
Bvec[x_,y_,z_]={Bvecx[x,y,z],Bvecy[x,y,z],Bvecz[x,y,z]};
Simplify[
Div[Avec[x,y,z]\[Cross]Bvec[x,y,z]]==Rot[Avec[x,y,z]].Bvec[x,y,z]-Avec[x,y,z].Rot[Bvec[x,y,z]]
]

Potential and its gradient

Force and electric field are negative gradients of the potential energy and electric potential, respectively:

F≡-∇U,         E≡-∇ϕ.

Gradient is perpendicular to the equipotential lines and it shows the direction of the strongest increase of the potential. (Thus the force shows the direction of the strongest decrease of the potential. This is natural because all systems tend to decrease their energy and the forces are developed accordingly to this). To prove this, consider an equipotential surface

U[x,y,z]=(Subscript[U, 0])

that goes through the point (Subscript[r, 0])={(Subscript[x, 0]), (Subscript[y, 0]),(Subscript[z, 0])}. Expanding the potential around (Subscript[r, 0]) up to the first order, one obtains

U[x,y,z]=(Subscript[U, 0])+(Subscript[U, x]^′)(x-(Subscript[x, 0]))+(Subscript[U, y]^′)(y-(Subscript[y, 0]))+(Subscript[U, z]^′)(z-(Subscript[z, 0])).

Together with the preceding formula this yields the equation of the plane tangential to the equipotential surface at (Subscript[r, 0])

(Subscript[U, x]^′)(x-(Subscript[x, 0]))+(Subscript[U, y]^′)(y-(Subscript[y, 0]))+(Subscript[U, z]^′)(z-(Subscript[z, 0]))=0,

where (Subscript[V, α]^′)≡(∂_α)V. Any vector r within this plane can be represented in the form

r=(Subscript[r, 0])+ρ,

where

ρ=(Subscript[e, x])(x-(Subscript[x, 0]))+(Subscript[e, y])(y-(Subscript[y, 0]))+(Subscript[e, z])(z-(Subscript[z, 0])).

One can see that  ρ is perpendicular to the gradient because

ρ·∇U=((Subscript[e, x])(x-(Subscript[x, 0]))+(Subscript[e, y])(y-(Subscript[y, 0]))+(Subscript[e, z])(z-(Subscript[z, 0])))·((Subscript[e, x])(Subscript[U, x]^′)+(Subscript[e, y])(Subscript[U, y]^′)+(Subscript[e, z])(Subscript[U, z]^′))
=(Subscript[U, x]^′)(x-(Subscript[x, 0]))+(Subscript[U, y]^′)(y-(Subscript[y, 0]))+(Subscript[U, z]^′)(z-(Subscript[z, 0]))=0.

As an example consider the electric potential created by two charges put at (0,0,0) and (a,0,0). In the (x,y) plane (z=0) the potential has the form

(Subscript[Q, 1])=1;   (Subscript[Q, 2])=-1;     a=1;
V[x_,y_]=((Subscript[Q, 1])/((((x^2)+(y^2)))^(1/2)))+((Subscript[Q, 2])/((((((x-a))^2)+(y^2)))^(1/2)));
Plot3D[V[x,y],{x,-2,2+a},{y,-2,2},PlotRange->{-5,5}]

EquipotentialLines=ContourPlot[V[x,y],{x,-2,2+a},{y,-2.5,2.5},Contours->20]

The electric field in the (x,y) plane is given by.

EE[x_,y_]=-Grad[V[x,y]]

Lines of the electric field are shown in the StreamPlot where we plot only x- and y-components

LinesE=StreamPlot[{EE[x,y][[1]],EE[x,y][[2]]},{x,-2,2+a},{y,-2.5,2.5},VectorPoints->10];

Combine the two plots

Show[EquipotentialLines,LinesE]

Lines of electric field should be perpendicular to the equipotential lines because they show the direction of the strongest decrease of the electric potential. However, it seems that Mathematica is not sufficiently precise in plotting to see it.