pythonoptimizationscipyfsolve

fsolve solver not converging to desired solutions


I try to find points on the surface of an ellipsoid that are also within a conical frustum. The ellipsoid and frustum are centred at ce and cf, respectively. The ellipsoid has semi-axes lengths defined in l. The frustum had a radius r1 at its base, an opening div, and a height h. To solve the issue at hand, I use python, together with numpy and scipy. More specifically, I use the fsolve function from scipy.optimize to minimise the difference between the ellipsoid and frustum equations. Based on that function, I have a Jacobian. My code is

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()

When the code is run, the points to which the code converges are the initial condition [t, p, z_val] or, if I try to put x0 as initial guess, the code converges to that value (which is not a valid solution), repeated over and over again. The solutions should be all the points of the frustum lying on the ellipsoid surface. I do not get why my code does not behave as expected? Can someone help me find out?


Solution

  • There are a number of problems with your approach.

    1. The intersection of the two surfaces is a curve. A curve in space has just one parameter, so you should not be using a triple-nested loop to do anything.

    2. You appear to be under the impression that the polar angle theta in the ellipsoid and the frustrum is the same. It is not. Theta and phi in the ellipsoid are just parameters: they are not actually the spherical polar angles (unless the ellipsoid happens to reduce to a sphere).

    3. There may be no intersections or there may in fact be two intersection curves, depending on the height of the frustrum. fsolve would never be able to find both (or even indicate which it had found).

    4. You can solve the problem exactly. For each polar angle phi for the generators of the frustrum there may be 0, 1 or 2 values of z, according to the number of intersections. These values are simply the real roots of a quadratic equation, which may be set up for each value of phi - see the bottom of this post.

    I give a code solution below and a sketch proof of the setting up of the quadratic equations at the bottom of this post. Note that I do not know whether you are referring to r1 as the larger or smaller radius of the frustrum. I have taken it as the larger one, but had to increase it to get any intersection. I do not use fsolve at all: there is no need and it could not solve the problem.

    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()
    

    Output: enter image description here

    Sketch summary of the solution. You will have to complete the lengthy algebra to find the coefficients of the quadratic: they are given in the code.

    enter image description here