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

Question

我需要计算网格上的值定义的表面下两点 (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()

标签： pythonnumpyscipyintegrationcalculus

Answer 1

如果你有表面点（我们甚至可以放宽规则网格的要求）和曲线点，那么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)

完整的设置如下所示：

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)

然后我们创建一个多维线性插值器：

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)

我们还计算二维曲线的曲线坐标：

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

使用基于累积弧长的曲线坐标，求解的积分如下：

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

最后我们使用我们选择的方法进行积分，这里是最简单的梯形规则：numpy

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