Line integral under surface defined by meshgrid values - Python

Question

I need to calculate the line integral between two points (x1,y1) and (x2,y2) under a surface defined by values on a meshgrid.

I'm not exactly sure on the best tool/approach to use for this process using python.

As I do not have a function which represents the surface, instead values at points on a evenly spaaced meshgrid I am assuming I will need to use one of the following methods

   trapz         -- Use trapezoidal rule to compute integral from samples.
   cumtrapz      -- Use trapezoidal rule to cumulatively compute integral.
   simps         -- Use Simpson's rule to compute integral from samples.
   romb          -- Use Romberg Integration to compute integral from
                    (2**k + 1) evenly-spaced samples.

Any help or guidance would be appreciated.

Edit:

import numpy as np
from scipy import interpolate

def f(x, y):
    return x**2 + x*y + y*2 + 1

xl = np.linspace(-1.5, 1.5, 101,endpoint = True)
X, Y = np.meshgrid(xl, xl)
Z = f(X, Y)

#And a 2D Line:
arr_2D = np.linspace(start=[-1, 1.2], stop=[0, 1.5], num=101,endpoint = 
True) #Creates a 2D line between these two points

#Then we create a multidimensional linear interpolator:
XY = np.stack([X.ravel(), Y.ravel()]).T
S = interpolate.LinearNDInterpolator(XY, Z.ravel())
print(S)

#To interpolate points from 2D curve on the 3D surface:
St = S(arr_2D)

#We also compute the curvilinear coordinates of the 2D curve:

#Using curvilinear coordinates based on cumulative arc length, the integral to solve looks like:
Sd = np.cumsum(np.sqrt(np.sum(np.diff(arr_2D, axis=0)**2, axis=1)))
print(Sd)

I = np.trapz(St[:-1], Sd) # 2.041770932394164
print("Integral: ",I)

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure()
ax = plt.axes(projection="3d")

x_line = np.linspace(start=[-1], stop=[1.5], num=100,endpoint = True)
y_line = np.linspace(start=[-1.2], stop=[1.5], num=100,endpoint = True)

ax.plot3D(x_line, y_line, 'red')  #Line which represents integral
ax.plot_wireframe(X, Y, Z, color='green') #Represents the surface
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Time')

plt.show()

fig = plt.figure()
ax = plt.axes()
ax.fill_between(Sd, St)
ax.set_xlabel('x')
ax.set_ylabel('Z')
plt.show()

Answer 1

Provided you have surface points (we can even relax the requirement of regular grid) and curve points, then basic analysis provided by numpy and scipy packages should do the trick.

First, let's create a trial dataset for your problem.

import numpy as np
from scipy import interpolate

Mainly a 3D surface:

def f(x, y):
    return x**2 + x*y + y*2 + 1

xl = np.linspace(-1.5, 1.5, 101)
X, Y = np.meshgrid(xl, xl)
Z = f(X, Y)

And a 2D curve:

t = np.linspace(0, 1, 1001)
xt = t**2*np.cos(2*np.pi*t**2)
yt = t**3*np.sin(2*np.pi*t**3)

The complete setup looks like:

axe = plt.axes(projection='3d')
axe.plot_surface(X, Y, Z, cmap='jet', alpha=0.5)
axe.plot(xt, yt, 0)
axe.plot(xt, yt, St)
axe.view_init(elev=25, azim=-45)

Then we create a multidimensional linear interpolator:

XY = np.stack([X.ravel(), Y.ravel()]).T
S = interpolate.LinearNDInterpolator(XY, Z.ravel())

To interpolate points from 2D curve on the 3D surface:

xyt = np.stack([xt, yt]).T
St = S(xyt)

We also compute the curvilinear coordinates of the 2D curve:

Sd = np.cumsum(np.sqrt(np.sum(np.diff(xyt, axis=0)**2, axis=1)))

Using curvilinear coordinates based on cumulative arc length, the integral to solve looks like:

fig, axe = plt.subplots()
axe.plot(Sd, St[:-1])
axe.fill_between(Sd, St[:-1], alpha=0.5)
axe.grid()

Finally we integrate using the method of our choice, here the simplest Trapezoidal Rule from numpy:

I = np.trapz(St[:-1], Sd) # 2.041770932394164