由网格值定义的曲面下的线积分--Python

问题描述

我需要计算由网格上的值定义的曲面下两点(x1，y1)和(x2，y2)之间的线积分。

我不太确定使用python进行此过程所使用的最佳工具/方法。

因为我没有表示曲面的函数，而是在均匀分布的网格上的点上的值，所以我假设我将需要使用以下方法之一

   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.

如有任何帮助或指导，我们将不胜感激。

编辑：

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()

---

解决方案

只要您有曲面点(我们甚至可以放宽规则网格的要求)和曲线点，则numpy和scipy包提供的基本分析应该可以完成此任务。

首先，让我们为您的问题创建一个试验数据集。

import numpy as np
from scipy import interpolate

主要是3D曲面：

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)

和二维曲线：

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)

完整设置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)

然后我们创建一个multidimensional linear interpolator：

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

从3D曲面上的2D曲线插值点：

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

我们还计算2D曲线的曲线坐标：

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

用curvilinear coordinates累加arc length求积分如下：

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

最后我们使用method of our choice进行集成，这里是numpy中最简单的Trapezoidal Rule：

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