fsolve求解不收敛：椭球与圆台交线问题排查

核心错误点 #

• 目标函数逻辑完全颠倒：原代码将theta, phi, z作为固定参数传入args，但实际需要求解的正是这三个变量——让椭球上的(theta, phi)点与圆台上的(theta, z)点重合。原函数本质上在求解恒等式，自然收敛到初始值。
• Jacobian矩阵求导对象错误：原矩阵的导数计算基于错误的变量定义，没有对实际待求解的theta, phi, z求偏导，无法为迭代提供有效梯度引导。
• 圆台参数范围错误：原代码中z的取值范围不符合frustum函数定义，该函数中z是圆台中心的轴向偏移量，正确范围应为0到圆台高度h，而非从cf[-1]开始。
• 冗余的三重循环：遍历所有theta, phi, z组合完全无效，交线在极角theta方向连续，只需遍历theta并求解对应phi和z即可。

修正后的代码 #

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

def get_ellipsoid(c, l):
    Theta, Phi = np.meshgrid(np.linspace(0, 2*np.pi, 50), np.linspace(0, np.pi, 50))
    X_ell = c[0] + l[0] * np.cos(Theta) * np.cos(Phi)
    Y_ell = c[1] + l[1] * np.sin(Theta) * np.cos(Phi)
    Z_ell = c[2] + l[2] * np.sin(Phi)
    return X_ell, Y_ell, Z_ell

def get_conical_frustum(c: tuple[float, float, float],
                        r1: float,
                        alpha: float,
                        h: float,
                    ) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
    Theta, Z = np.meshgrid(np.linspace(0, 2*np.pi, 50), np.linspace(0, h, 50))
    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_point(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_point(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
    ])

# 修正目标函数：待求解变量为[phi, z]，theta作为固定参数
def func(x, theta, ce, cf, l, r1, div):
    phi, z = x
    return ellipsoid_point(theta, phi, ce, l) - frustum_point(theta, z, cf, r1, div)

# 修正Jacobian：对phi和z求偏导
def jacobian(x, theta, ce, cf, l, r1, div):
    phi, z = x
    d_ellipsoid_dphi = np.array([
        -l[0] * np.cos(theta) * np.sin(phi),
        -l[1] * np.sin(theta) * np.sin(phi),
        l[2] * np.cos(phi)
    ])
    d_frustum_dz = np.array([
        np.tan(div) * np.cos(theta),
        np.tan(div) * np.sin(theta),
        1.0
    ])
    return np.column_stack([d_ellipsoid_dphi, -d_frustum_dz])

# 参数定义
ce = np.array([0.0, 0.0, 0.0])
l = np.array([1.5, 1.1, 0.45])

cf = np.array([0.0, 0.0, -1.5])
h = np.linalg.norm(ce - cf)
r1 = 0.15
div = np.pi / 18

# 遍历theta求解交线点
num_theta = 50
theta_list = np.linspace(0, 2*np.pi, num_theta)
intersection_points = []

for theta_val in theta_list:
    x0 = np.array([0.0, h/2])  # 初始猜测：椭球赤道面+圆台中间高度
    try:
        sol = fsolve(func, x0, args=(theta_val, ce, cf, l, r1, div), fprime=jacobian)
        phi_sol, z_sol = sol
        if 0 <= z_sol <= h:  # 验证解在圆台范围内
            point = ellipsoid_point(theta_val, phi_sol, ce, l)
            intersection_points.append(point)
    except RuntimeError:
        continue

intersection_points = np.array(intersection_points)

# 绘图
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection='3d')

xell, yell, zell = get_ellipsoid(ce, l)
ax.plot_surface(xell, yell, zell, color='r', alpha=0.1)

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

if len(intersection_points) > 0:
    ax.scatter(intersection_points[:,0], intersection_points[:,1], intersection_points[:,2], color='g', s=10, label='Intersection')
ax.legend()

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

修正说明 #

• 目标函数重构：将待求解变量改为[phi, z]，theta作为固定参数，符合几何交线的单theta对应唯一解的逻辑。
• Jacobian修正：正确计算差值函数对phi和z的偏导数，为迭代提供有效梯度，引导收敛到真实解。
• 参数范围修正：调整圆台z的取值范围为0到h，匹配函数定义。
• 优化遍历逻辑：仅遍历theta，减少无效计算，同时增加解的有效性验证，过滤超出圆台范围的结果。

内容的提问来源于stack exchange，提问作者aheuchamps