Решатель fsolve не сходится к желаемым решениям

Ваш подход имеет ряд проблем.

1. Пересечение двух поверхностей представляет собой кривую. Кривая в пространстве имеет только один параметр, поэтому вам не следует использовать тройной вложенный цикл для каких-либо действий.

2. Похоже, у вас сложилось впечатление, что полярный угол тета в эллипсоиде и усеченной пирамиде один и тот же. Это не. Тета и фи в эллипсоиде — это просто параметры: на самом деле они не являются сферическими полярными углами (если только эллипсоид не превратится в сферу).

3. Пересечений может не быть или на самом деле могут быть две кривые пересечения, в зависимости от высоты усеченного конуса. fsolve никогда не сможет найти оба (или даже указать, что он нашел).

4. Вы можете решить проблему точно. Для каждого полярного угла фи для образующих пирамиды может быть 0, 1 или 2 значения z в зависимости от количества пересечений. Эти значения являются просто действительными корнями квадратного уравнения, которое можно задать для каждого значения фи — см. конец этого поста.

Ниже я привожу решение кода и эскизное доказательство построения квадратных уравнений внизу этого поста. Обратите внимание: я не знаю, имеете ли вы в виду 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()

Выход:

Краткое описание решения. Вам придется выполнить длинную алгебру, чтобы найти коэффициенты квадратичного: они заданы в коде.