2 Numerical Integrals

We can use the Taylor series to calculate numerical approximations to integrals too. Let's look at integrating the function f (x) from x₀ to x₀ + h numerically. Expanding f (x) around the midpoint of the interval we have

$$\int_{a}^{a+h} \int_{a}^{a+h} \left [ f(a+h/2) + f'(a+h/2) \left (x-a-\frac{h}{2} \right) +O(h^2)\right]d$$ =

= f (a + h/2)h + O(h ³) (13)

To apply this to the integral from a to b , we can divide it into N subintervals with

$${\frac{b-a}{N}}$$ h =

$$\left (i - \frac{1}{2} \right)~$$ x_i = (14)

and
 

$$\int_{a}^{b} \sum_{i=1}^{N}$$ =

$${\frac{b-a}{N}} \sum_{i=1}^{N}$$ = (15)

It is interesting to notice that round off will produce an error of approximately rf (x) in f (x) , where r is the precision defined in section 1, which bounded by (b - a)rmax (f (x)) , which is independent of N if the sum is done carefully. Integrations tend to be much more stable than derivatives when decreasing the step size.

Now that we know how to do a simple integral let's see how to do the various kinds of integrals we will see in electrodynamics. Often we need line integrals like

$$\int_{\vec \ell}^{} \vec{A} \cdot \vec{\ell}$$ (16)

either around a loop or along a path between two end points. One way is to first define the path parametrically. For example if we need to calculate the line integral of $\vec{A}$($\vec{r}$) around an ellipse with semimajor axis of length a parallel to $\hat{x}$ , semiminor axis of length b in the y - z plane making an angle of $\alpha$ with the y axis, we can define the path parametrically as

$$\vec{\ell} \hat{x} \alpha \hat{y} \alpha \hat{z}$$ (17)

where u runs from 0 to 2$\pi$ . The line integral can be approximated by
 

$$\equiv {\frac{2 \pi}{N}}$$ u_i

$$\vec{r}_{i}^{} \equiv \vec{\ell} \left (\frac{u_i+u_{i-1}}{2} \right)$$

$$\vec{h}_{i}^{} \equiv \vec{\ell} \vec{\ell}$$

$$\int_{\vec \ell}^{} \vec{A} \cdot \vec{\ell} \sum_{i=1}^{N} \vec{A} \vec{r}_{i}^{} \cdot \vec{h}_{i}^{}$$ = (18)

Notice that this is written in a form exactly analogous to Eq. 15, and is an explicit example of what is shown in figure 1.3 of Jackson. The generalization to other paths should now be easy.

Next we need to consider two-dimensional integrals. One way to perform these is to apply the one-dimensional method for each dimension, or equivalently Taylor series expand in both x and y . To integrate over a small rectangle of sides h_a and h_b ,

$$\int_{a}^{a+h_a} \int_{b}^{b+h_b} \int_{a}^{a+h_a} \int_{b}^{b+h_b} \left [ f(a+h_a/2,b+h_b/2) \right.$$ =

$${\frac{h_a}{2}}$$ +

$$\left . f_y(a+h_a/2,b+h_b/2) (y-b-\frac{h_b}{2}) +O(h^2) \right]d$$ +

= f (a + h_a/2,b + h_b/2)h_ah_b + O(h ⁴) (19)

where O(h ⁴) etc. means that order of products of the factors h_a and h_b , and f_x and f_y are the partial derivatives of f with respect to x and y . Notice this is nothing but the value of the function at the center of the rectangle multiplied by the area of the rectangle. If we integrate over an area using N points total (that is N little rectangles or roughly N ^1/2 in each direction), the integration error will be O(N ^{- 1}) . Notice we can easily generalize this to other basic shapes by multiplying a value of the function at a central point times the area of the shape.

Going through the same mathematics for a three dimensional integral, we get the obvious generalization of the value of the function at the center of the 3-d rectangular prism times the volume of the prism,

$$\int_{a}^{a+h_a} \int_{b}^{b+h_b} \int_{c}^{c+h_c}$$ (20)

and if we integrate over a volume using N points total, the integration error will be O(N ^{- 2/3}) . For a general d dimensional integration, the midpoint rule will give an error that goes like O(N ^{- 2/d}) . Again we can generalize to other shapes by multiplying a value of the function at a central point times the volume of the shape.

It is instructive to do the line integral of $\vec{E}$ around the outside of a very small square of side h centered at (x₀,y₀,z₀) . Let's take the square to be oriented along the x and y axes. The line integral is,

$$\int \vec{E} \vec{r} \cdot \vec{\ell}$$ (21)

Notice that the right hand side is our central difference formula for the z component of the curl at the center of the square, multiplied by h ². That is, we have proved Stokes theorem. A similar application to a surface integral over a little cube will give the divergence theorem.

Let's see how to apply our integration to a surface integral like those typically found in electromagnetism. Often we want to integrate the surface normal component of a vector field over a surface. For example, we might want to integrate the normal component of the electric field over a surface

$$\int \vec{E} \cdot \vec{A}$$ (22)

Rather than do a trivial problem like the surface of a cube or sphere, lets integrate over an ellipsoidal surface. We do this problem in two ways. First we imagine we are going to set up the problem analytically and then do the double integral numerically at the end. Second we will imagine doing it numerically from the beginning.

One way to define the surface is parametrically, similarly to the way we defined a path parametrically. We pick cartesian coordinates so that the origin is at the center of our ellipsoid with the coordinate axes aligned with the semiaxes of the ellipsoid. The surface can be given by

x(u,v) = asin(u)cos (v)

y(u,v) = bsin(u)sin (v)

z(u,v) = ccos (u) (23)

with 0 $\leq$ u $\leq$ $\pi$ , 0 $\leq$ v $\leq$ 2$\pi$ . As with any parametric definition, there are an infinity of choices you can make. This parameterization defines the vector

$$\vec{S} \hat{x} \hat{y} \hat{z}$$ (24)

which describes our surface. The directed area of the surface corresponding to a small increment in u and v is

$$\vec{A} \left [ \frac{\partial \vec S(u,v)}{\partial u} \times \frac{\partial \vec S(u,v)}{\partial v} \right]d$$ (25)

where we have chosen the sign for our particular problem so that d$\vec{A}$ points outward. The analytic expression for the surface integral is

$$\int_{0}^{2\pi} \int_{0}^{\pi} \vec{E} \cdot \left [ \frac{\partial \vec S(u,v)}{\partial u} \times \frac{\partial \vec S(u,v)}{\partial v} \right]$$ (26)

Note: if you have trouble following this to here, you are having trouble with standard analytic vector calculus not numerical vector calculus. To attempt an analytic integration, you would do the derivatives above and simplify the expression.

We will now convert our expression into a form ready for numerical calculation. We could do the derivatives analytically, but for amusement lets do the integrals using the midpoint rule and the derivatives using the central differences.
 

$${\frac{\pi}{N_u}}$$ h_u =

$${\frac{2 \pi}{N_v}}$$ h_v =

$${\textstyle\frac{1}{2}}$$ u_i =

$${\textstyle\frac{1}{2}}$$ v_j =

$$\Delta \vec{A}_{ij}^{} \left [ \vec S(u_i+h_u/2,v_j) - \vec S(u_i-h_u/2,v_j) \right]\times \left [ \vec S(u_i,v_j+h_v/2) - \vec S(u_i,v_j-h_v/2) \right]$$ =

$$\int \vec{E} \cdot \vec{A} \sum_{i=1}^{N_u} \sum_{j=1}^{N_v} \vec{E} \cdot \Delta \vec{A}_{ij}^{}$$ = (27)

Alternatively, you can think about doing this surface integral from a numerical point of view from the beginning. In that case, we could break up u and v into equally spaced intervals. This would divide the surface into little 4 sided areas (actually at u = 0 and u = $\pi$ one of the sides goes to zero and they become 3 sided there, but this can be ignored for now) and we can think of doing the integration by multiplying the surface area made up of the little areas bounded by the four surface points

$$\vec{S}_{i,j,\pm,\pm}^{} \vec{S} \pm \pm$$ (28)

If we approximate this area by that of a straight sided quadrilateral with an outward direction we get the cross products of the little vectors that go across the quadrilateral or an equivalent expression. This is just the cross product given in Eq. 27 above. Notice that we can use this ``numerical'' thinking to derive or check our analytic results and gain a better understanding of how everything fits together.