由网格值定义的曲面下的线积分 - Python
Posted
技术标签:
【中文标题】由网格值定义的曲面下的线积分 - Python【英文标题】:Line integral under surface defined by meshgrid values - Python 【发布时间】:2020-12-25 08:14:03 【问题描述】:我需要计算网格网格上的值定义的曲面下两点 (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()
【问题讨论】:
【参考方案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)
完整设置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)
我们还计算二维曲线的曲线坐标:
Sd = np.cumsum(np.sqrt(np.sum(np.diff(xyt, axis=0)**2, axis=1)))
使用基于累积arc length 的curvilinear coordinates,求解的积分如下:
fig, axe = plt.subplots()
axe.plot(Sd, St[:-1])
axe.fill_between(Sd, St[:-1], alpha=0.5)
axe.grid()
最后我们使用method of our choice进行集成,这里最简单的Trapezoidal Rule来自numpy
:
I = np.trapz(St[:-1], Sd) # 2.041770932394164
【讨论】:
谢谢你帮了我很多忙!我正在尝试使用曲线坐标重新创建您在底部生成的线积分图。我究竟做错了什么 ? fig = plt.figure() ax = plt.axes() ax.fill_between(Sd, St) ax.set_xlabel('x') ax.set_ylabel('Z') plt.show() 这很有帮助。如果您有兴趣,我仍在努力解决此处解释的问题link @Morgan,好问题,现在我明白你为什么需要解决这个小问题了。 @jilandercy 可以修改上述函数以计算投影平行六面体上的积分吗? @jlandercy 是的,我想我没有意识到您在这里使用的是函数而不是真实数据。我已经设法使用 griddata、RectBavariateSpline 和积分函数解决了我的问题 :)以上是关于由网格值定义的曲面下的线积分 - Python的主要内容,如果未能解决你的问题,请参考以下文章