fsolve 求解器未收敛到所需的解

Question

我尝试在椭球体表面上找到位于截头圆锥体内的点。椭球体和平截头体分别以

ce

和

cf

为中心。椭球体的半轴长度在

中定义。平截头体的底部具有半径

r1

、开口

div

和高度

。为了解决手头的问题，我使用

python

以及

numpy

和

scipy

。更具体地说，我使用

fsolve

中的

scipy.optimize

函数来最小化椭球体方程和平截头体方程之间的差异。基于该函数，我有一个雅可比行列式。我的代码是

import numpy as np
from scipy.optimize import fsolve
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

def get_ellispoid(c, l):
    X_ell = c[0] + l[0] * np.outer(np.sin(Phi), np.cos(Theta))
    Y_ell = c[1] + l[1] * np.outer(np.sin(Phi), np.sin(Theta))
    Z_ell = c[2] + l[2] * np.outer(np.cos(Phi), np.ones_like(Theta))
    return X_ell, Y_ell, Z_ell

def get_conical_frustum(c: tuple[float, float, float],
                        r1: float,
                        alpha: float,
                        z: np.ndarray,
                    ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    Theta, Z = np.meshgrid(theta, z)
    R = r1 + Z * np.tan(alpha)
    X_fru = c[0] + R * np.cos(Theta)
    Y_fru = c[1] + R * np.sin(Theta)
    Z_fru = c[2] + Z
    return X_fru, Y_fru, Z_fru

def ellipsoid(theta, phi, ce, l):
    return np.array(
                [ce[0] + l[0] * np.cos(theta) * np.cos(phi),
                    ce[1] + l[1] * np.sin(theta) * np.cos(phi),
                    ce[2] + l[2] * np.sin(phi)]
            )

def frustum(theta, z, cf, r1, div):
    return np.array(
                [cf[0] + (r1 + z * np.tan(div)) * np.cos(theta),
                    cf[1] + (r1 + z * np.tan(div)) * np.sin(theta),
                    cf[2] + z]
            )


# Define the function for which we want to find the roots
def func(x, theta, phi, z, ce, cf, l, r1, div):
    return ellipsoid(theta, phi, ce, l) - frustum(theta, z, cf, r1, div)

# Define the Jacobian matrix of the function
def jacobian(x, theta, phi, z, ce, cf, l, r1, div):
    return np.array(
        [
            [-(l[0] * np.cos(phi) - r1 - z * np.tan(div)) * np.sin(theta), -l[0] * np.cos(theta) * np.sin(phi), -np.tan(div) * np.cos(theta)],
            [-(l[1] * np.cos(phi) - r1 - z * np.tan(div)) * np.cos(theta), -l[1] * np.cos(theta) * np.sin(phi), -np.tan(div) * np.sin(theta)],
            [0.0, l[2] * np.cos(phi), -1],
        ]
    )

num_points = 20
theta = np.linspace(0, 2*np.pi, num_points)
phi = np.linspace(0, np.pi, num_points)
Theta, Phi = np.meshgrid(theta, phi)


# Ellipsoid
ce = [0.0, 0.0, 0.0]
l = [1.5, 1.1, 0.45]

# Conical frustum
cf = [0.0, 0.0, -1.5]
h = np.linalg.norm(np.array(cf) - np.array(ce))
z = np.linspace(cf[-1], h, num_points)
r1 = 0.15
div = np.pi / 18

# Initial guess
x0 = np.array([0.0, 0.0, ce[-1]])

# Iterate over theta, phi, and z values
sol = []
for it, t in enumerate(theta):
    for ip, p in enumerate(phi):
        for iz, z_val in enumerate(z):
            try:
                root = fsolve(func, [t, p, z_val], args=(t, p, z_val, ce, cf, l, r1, div), fprime=jacobian)
                sol.append(root)
            except RuntimeError:
                continue
sol = np.array(sol)

# Create a 3D plot
fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection = '3d')

# Plot the volumes
xell, yell, zell = get_ellispoid(ce, l)
ax.plot_surface(xell, yell, zell, color = 'r', alpha = 0.005)

xfru, yfru, zfru = get_conical_frustum(cf, r1, div, z)
ax.plot_surface(xfru, yfru, zfru, color = 'b', alpha = 0.35)

# Plot the intersections
ax.scatter(sol[:, 0], sol[:, 1], sol[:, 2], color = 'g', s = 1.5)

# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

# Show the plot
plt.show()

当代码运行时，代码收敛的点是初始条件

[t, p, z_val]

，或者，如果我尝试将

x0

作为初始猜测，代码收敛到该值（这不是有效的解决方案），一遍又一遍地重复。解应该是位于椭球面上的平截头体的所有点。我不明白为什么我的代码没有按预期运行？有人可以帮我查一下吗？

Answer 1

您的方法存在很多问题。

两个曲面的交线是一条曲线。空间中的曲线只有一个参数，因此您不应该使用三重嵌套循环来执行任何操作。
您似乎认为椭球体和平截头体中的极角 θ 是相同的。它不是。椭球体中的 Theta 和 phi 只是参数：它们实际上并不是球极角（除非椭球体恰好缩小为球体）。
可能没有相交，或者实际上可能有两条相交曲线，具体取决于平截头体的高度。 fsolve 永远无法找到两者（甚至无法指出它找到了哪一个）。
你可以准确地解决问题。对于平截头体生成元的每个极角 phi，根据交点的数量，可能有 0、1 或 2 个 z 值。这些值只是二次方程的实根，可以为 phi 的每个值设置 - 请参阅本文的底部。

我在下面给出了一个代码解决方案，并在本文底部给出了二次方程设置的草图证明。请注意，我不知道您是否将 r1 称为截锥体的较大半径或较小半径。我已将其视为较大的一个，但必须增加它才能获得任何交叉点。我根本不使用 fsolve：没有必要，它无法解决问题。

import math
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# Ellipsoid
xe, ye, ze = 0.0, 0.0, 0.0
a, b, c = 1.5, 1.1, 0.45

# Conical frustum
xf, yf, zf = 0.0, 0.0, -1.5
h = 1.5
R = 0.45
alpha = math.pi / 18
tanalpha = math.tan( alpha )

angles = np.linspace( 0.0, 2 * np.pi, 40 )
upper = []
lower = []

for phi in angles:
   # Set up coefficients in quadratic equation
   cosphi, sinphi = math.cos( phi ), math.sin( phi )
   x1, x2 = xf - xe + ( R + zf * tanalpha ) * cosphi, tanalpha * cosphi
   y1, y2 = yf - ye + ( R + zf * tanalpha ) * sinphi, tanalpha * sinphi
   z1, z2 = ze, 1
   A = ( x2 / a ) ** 2 + ( y2 / b ) ** 2 + ( z2 / c ) ** 2
   B = -2 * ( x1 * x2 / a ** 2 + y1 * y2 / b ** 2 + z1 * z2 / c ** 2 )
   C = ( x1 / a ) ** 2 + ( y1 / a ) ** 2 + ( z1 / c ) ** 2 - 1
   discriminant = B ** 2 - 4 * A * C
   if discriminant >= 0:
      root1 = ( -B + math.sqrt( discriminant ) ) / ( 2 * A )
      root2 = -B / A - root1
      if zf <= root1 <= zf + h: lower.append( [ xf + ( R - ( root1 - zf ) * tanalpha ) * cosphi, yf + ( R - ( root1 - zf ) * tanalpha ) * sinphi, root1 ] )
      if zf <= root2 <= zf + h: upper.append( [ xf + ( R - ( root2 - zf ) * tanalpha ) * cosphi, yf + ( R - ( root2 - zf ) * tanalpha ) * sinphi, root2 ] )

lower = np.array( lower )
upper = np.array( upper )

################################################################################################

# Create a 3D plot

def get_ellipsoid( c, l ):
    num_points = 20
    th = np.linspace( 0, 2*np.pi, num_points )
    ph = np.linspace( 0, np.pi, num_points )
    Theta, Phi = np.meshgrid(th, ph)

    X_ell = c[0] + l[0] * np.outer(np.sin(Phi), np.cos(Theta))
    Y_ell = c[1] + l[1] * np.outer(np.sin(Phi), np.sin(Theta))
    Z_ell = c[2] + l[2] * np.outer(np.cos(Phi), np.ones_like(Theta))
    return X_ell, Y_ell, Z_ell

def get_frustum( c, r1, alpha, h ):
    num_points = 20
    th = np.linspace( 0, 2*np.pi, num_points )
    z = np.linspace( 0, h, num_points )
    Theta, Z = np.meshgrid(th, z)
    R = r1 - Z * np.tan(alpha)
    X_fru = c[0] + R * np.cos(Theta)
    Y_fru = c[1] + R * np.sin(Theta)
    Z_fru = c[2] + Z
    return X_fru, Y_fru, Z_fru




fig = plt.figure(figsize = (8, 6))
ax = fig.add_subplot(111, projection = '3d')

# Plot the volumes
xell, yell, zell = get_ellipsoid( [xe,ye,ze], [a,b,c] )
ax.plot_surface(xell, yell, zell, color = 'r', alpha = 0.005)

xfru, yfru, zfru = get_frustum( [xf,yf,zf], R, alpha, h )
ax.plot_surface( xfru, yfru, zfru, color = 'b', alpha = 0.35)

# Plot the intersections
if lower.size > 0: ax.plot( lower[:,0], lower[:, 1], lower[:, 2], color = 'g' )
if upper.size > 0: ax.plot( upper[:,0], upper[:, 1], upper[:, 2], color = 'g' )

# Set labels and title
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')

plt.show()

输出：

解决方案的概要。您必须完成冗长的代数才能找到二次方程的系数：它们在代码中给出。

fsolve 求解器未收敛到所需的解

问题描述投票：0回答：1

1个回答

最新问题

fsolve 求解器未收敛到所需的解

问题描述 投票：0回答：1

1个回答

最新问题

问题描述投票：0回答：1