Example 1
A line integral is an integral of a function along a curve in
space. We usually represent the curve by a parametric equation,
e.g.
Then, we can write the line integral definition as:
where
The following examples are adapted from Chapter 10 in Advanced Engineering Mathematics by Kreysig.
The first example is the evaluation of a line integral in the
plane. We want to evaluate the integral of
import autograd.numpy as np
from autograd import jacobian
from scipy.integrate import quad
def F(X):
x, y = X
return -y, -x * y
def r(t):
return np.array([-np.sin(t), np.cos(t)])
drdt = jacobian(r)
def integrand(t):
return F(r(t)) @ drdt(t)
I, e = quad(integrand, 0.0, np.pi / 2)
print(f'The integral is {I:1.4f}.')
The integral is 0.4521.
We get the same result as the analytical solution.
The next example is in three dimensions. Find the line integral along
import autograd.numpy as np
from autograd import elementwise_grad, grad, jacobian
def F(X):
x, y, z = X
return [z, x, y]
def C(t):
return np.array([np.cos(t), np.sin(t), 3 * t])
dCdt = jacobian(C, 0)
def integrand(t):
return F(C(t)) @ dCdt(t)
I, e = quad(integrand, 0, 2 * np.pi)
print(f'The integral is {I:1.2f}.')
The integral is 21.99.
That is also the same as the analytical solution. Note the real analytical solution was 7 π, which is nearly equivalent to our answer.
print (7 * np.pi - I)
3.552713678800501e-15
As a final example, we consider an alternate form of the line integral. In this form we do not use a dot product, so the integral results in a vector. This doesn't require anything from autograd, but does require us to be somewhat clever in how to do the integrals since quad can only integrate scalar functions. We need to integrate each component of the integrand independently. Here is one approach where we use lambda functions for each component. You could also manually separate the components.
def F(r):
x, y, z = r
return x * y, y * z, z
def r(t):
return np.array([np.cos(t), np.sin(t), 3 * t])
def integrand(t):
return F(r(t))
[quad(lambda t: integrand(t)[i], 0, 2 * np.pi)[0] for i in [0, 1, 2]]
Out[1]: [-6.845005548807377e-17, -18.849555921538755, 59.21762640653615]
[0, -6 * np.pi, 6 * np.pi**2]
Out[1]: [0, -18.84955592153876, 59.21762640653615]
which is evidently the same as our numerical solution.
Maybe an alternative, but more verbose is this vectorized integrate function. We still make temporary functions for integrating, and the vectorization is essentially like the list comprehension above, but we avoid the lambda functions.
@np.vectorize
def integrate(i):
def integrand(t):
return F(r(t))[i]
I, e = quad(integrand, 0, 2 * np.pi)
return I
integrate([0, 1, 2])
Out[1]: array([-6.84500555e-17, -1.88495559e+01, 5.92176264e+01])