In this example we illustrate Gauss's theorem, Green's identities, and Stokes' theorem in Chebfun3.
1. Gauss's theorem
Gauss's theorem, also known as the divergence theorem, asserts that the integral of the sources of a vector field in a domain $K$ is equal to the flux of the vector field through the boundary: $$ \int_K \mbox{div}(\vec{v}) \, dV = \int_{\partial K} \vec{v} \cdot \vec{dS}. $$ Here $\vec{dS}$ is the vectorial surface element given by $\vec{dS} = \vec{n} dS$, where $\vec{n}$ is the outward normal vector to the surface $\partial K$ and $dS$ is the surface element.
Let us consider a simple 3D vector field $\vec{v}$ on the cube $K = [-1, 1]^3$:
cheb.xyz; v1 = x.^2-y; v2 = y.^2; v3 = z; v = [ v1 ; v2 ; v3 ];
Here we have built $\vec{v}$ by stacking three Chebfun3 objects one over another, beginning with the command cheb.xyz
to generate chebfun3 objects x
, y
, and z
corresponding to the functions $x$, $y$, and $z$. Alternatively we could pass three function handles directly to the Chebfun3v constructor; see the example for Stokes' Theorem below.
The left-hand side of the identity of Gauss's theorem, the integral of the divergence, can be computed in Chebfun3 like this, nicely matching the exact value 8:
format long I1 = sum3(div(v))
I1 = 7.999999999999998
For the right-hand side we need the flux integral over the boundary of the cube, which consists of six surfaces. The outward normal vectors on the sides ${ \pm 1 } \times [-1, 1]^2$ are $\pm [1; 0; 0]$. Thus we have to integrate $\vec{v}(\pm 1, y, z) \cdot [\pm 1; 0; 0] = \pm v_1(\pm 1, y, z)$ over these sides of the cube. The restrictions $v_1(\pm 1,:,:)$ are represented as chebfun2 objects. Treating the other two pairs of sides similarly, we compute a result again with excellent accuracy:
I2 = sum2(v1(1,:,:)) - sum2(v1(-1,:,:)) + ... sum2(v2(:,1,:)) - sum2(v2(:,-1,:)) + ... sum2(v3(:,:,1)) - sum2(v3(:,:,-1))
I2 = 8
2. Green's identities
The first Green identity is an analogue of integration by parts in higher dimensions: $$ \int_K ( f \Delta g + \nabla f \cdot \nabla g ) \, dV = \int_{\partial K} f \nabla g \cdot \vec{dS}. $$ This is a corollary of the Gauss's theorem (set $\vec{v} = f \nabla g$). The second Green identity is $$ \int_K (f \Delta g - (\Delta f) g ) \, dV = \int_{\partial K} ( f \nabla g - g \nabla f ) \cdot \vec{dS}. $$ In both formulas $f$ and $g$ are scalar functions on $K$. To illustrate both formulas in Chebfun3, let us consider the functions
f = 1 + x.*exp(y+z); g = x.^2 + y.^2 + z.^2; isosurface(g,1.3,'r'), axis equal, grid on
The Laplacian and gradient of $f$ can be computed with lap(f)
and grad(f)
. We evaluate the two sides of the first Green identity in the same way as for the Gauss's theorem above. The exact value $48$ is matched closely by the integral over the cube,
I3 = sum3(f .* lap(g) + dot(grad(f), grad(g)))
I3 = 48.000000000000114
The flux integral also agrees:
v = f * grad(g); v1 = v(1); v2 = v(2); v3 = v(3); I4 = sum2(v1(1,:,:)) - sum2(v1(-1,:,:)) + ... sum2(v2(:,1,:)) - sum2(v2(:,-1,:)) + ... sum2(v3(:,:,1)) - sum2(v3(:,:,-1))
I4 = 47.999999999999993
For the second Green formula the exact value of the integrals is again $48$. We compute
I5 = sum3(f .* lap(g) - lap(f) .* g)
I5 = 47.999999999999773
and
v = f * grad(g) - g * grad(f); v1 = v(1); v2 = v(2); v3 = v(3); I6 = sum2(v1(1,:,:)) - sum2(v1(-1,:,:)) + ... sum2(v2(:,1,:)) - sum2(v2(:,-1,:)) + ... sum2(v3(:,:,1)) - sum2(v3(:,:,-1))
I6 = 47.999999999999886
3. Stokes' Theorem
Stokes' Theorem in its classical formulation takes the form $$ \int_S \mbox{curl}(\vec{v}) \cdot \vec{dS} = \int_{\partial S} \vec{v} \cdot \vec{ds}. $$ Here, the vectorial line element $\vec{ds}$ is $\vec{ds} = \vec{t} ds$, where $\vec{t}$ is the tangential vector and $ds$ is the scalar line element. As an easy example we consider the surface given by the unit disk in the $x$ - $y$ plane, which can be parametrised by
S = chebfun2v(@(rho,phi) rho.*cos(phi), @(rho,phi) rho.*sin(phi), ... @(rho,phi) 0, [0, 1, 0, 2*pi]); surf(S)
We consider again the vector field $\vec{v}$ from the beginning, which we now construct by passing three function handles to the Chebfun3v constructor:
v = chebfun3v(@(x,y,z) x.^2 - y, @(x,y,z) y.^2, @(x,y,z) z); quiver3(v)
The curl of $\vec{v}$ is computed by
curlv = curl(v); quiver3(curlv)
Note how the curl of $\vec{v}$ points up, just through the surface we are considering. The flux integral then is
I7 = integral2(curlv, S)
I7 = 3.141592653589792
The exact value is $\pi$,
pi
ans = 3.141592653589793
To compute the line integral over the boundary of the surface, we first parametrize the circle by
gamma = chebfun(@(t) [cos(t), sin(t), 0*t], [0, 2*pi]);
Then the line integral is
I8 = integral(v, gamma)
I8 = 3.141592653589793