使用SciPy最小化解决悬链线问题

Question

我正在尝试使用 SciPy 优化来找到最小化两点之间悬挂绳索势能的函数，即悬链线问题。目标函数定义为

由于绳索的长度是固定的，因此弧长必须等于某个固定值。

我尝试使用以下命令将其编程到 SciPy 中：

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import minimize

y0 = 1
# Example: Define the x domain and initial guess for y
x_vals = np.linspace(0, 10, 100)  # x values from 1 to 10, 100 points
initial_y = np.ones_like(x_vals)  # Initial guess for y(x)
def objective(y_vals):
    # Number of points
    n = len(x_vals)
    
    # Reshape y_vals to be the correct shape for our problem
    y = y_vals
    
    # Calculate the finite differences for dy/dx
    dy_dx = np.diff(y) / np.diff(x_vals)
    
    # The objective function is the integral, which we approximate with the sum of the terms
    integrand = y_vals[:-1] * np.sqrt(1 + dy_dx**2)  # Use y[:-1] since dy_dx has one less point
    
    # Discretize the integral: sum the contributions from each segment
    integral = np.sum(integrand * (x_vals[1:] - x_vals[:-1]))
    
    return integral

def arc(y_vals):
    # Number of points
    n = len(x_vals)
    
    # Reshape y_vals to be the correct shape for our problem
    y = y_vals
    
    # Calculate the finite differences for dy/dx
    dy_dx = np.diff(y) / np.diff(x_vals)
    
    integrand = np.sqrt(1 + dy_dx**2)
    integral = np.sum(integrand * (x_vals[1:] - x_vals[:-1]))
    
    return integral

# Boundary conditions
def boundary_conditions(y_vals):
    return [y_vals[0] - y0, y_vals[-1]-y0, arc(y_vals)-12]  # set y0 at start and end and arc length



# Define the constraints for the optimization
constraints = [{'type': 'eq', 'fun': lambda y: boundary_conditions(y)}]

# Minimize the objective function with the constraints
result = minimize(objective, initial_y, constraints=constraints)

# Output the optimized values of y and the minimized objective value
optimized_y = result.x
print("Optimized y values:", optimized_y)
print("Minimized objective value:", result.fun)

plt.plot(optimized_y)
plt.show()

然而，结果却与悬链线原本应有的样子完全不同：

我怀疑这可能与我离散积分的方式有关，但我不确定。解决这个问题的正确方法是什么？

Answer 1

以下工作直至分辨率为 49：

import functools

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize, NonlinearConstraint, Bounds


def unpack(x: np.ndarray, ymid: np.ndarray, yend: float) -> tuple[
    np.ndarray,  # y including endpoints
    np.ndarray,  # x discrete differential
    np.ndarray,  # discrete derivative
    np.ndarray,  # integrand
]:
    y = np.empty(shape=ymid.size + 2, dtype=ymid.dtype)
    y[[0, -1]] = yend
    y[1: -1] = ymid
    dx = np.diff(x)
    dy = np.diff(y)
    dy_dx = dy/dx
    integrand = np.sqrt(1 + dy_dx**2)
    return y, dx, dy_dx, integrand


def potential_energy(ymid: np.ndarray, x: np.ndarray, yend: float) -> float:
    y, dx, dy_dx, integrand = unpack(x, ymid, yend)
    integrand *= y[:-1]  # Use y[:-1] since dy_dx has one less point
    return integrand.dot(dx)


def arc(ymid: np.ndarray, x: np.ndarray, yend: float) -> float:
    y, dx, dy_dx, integrand = unpack(x, ymid, yend)
    return integrand.dot(dx)


def solve(
    width: float = 10.,
    resolution: int = 101,
    length: float = 12.,
) -> tuple[np.ndarray, np.ndarray]:
    x = np.linspace(start=-0.5*width, stop=0.5*width, num=resolution)

    # The lower bound height is seen where the rope is pulled to two line segments forming
    # two right triangles; each has a hypotenuse of 0.5*length and a width of 0.5*width.
    lower_bound = -0.5*np.sqrt(length**2 - width**2)
    y0 = np.empty_like(x)
    y0[:1+resolution//2] = x[:1+resolution//2]*(2*lower_bound/width) + lower_bound
    y0[1+resolution//2:] = y0[resolution//2-1::-1]
    y0 -= lower_bound  # positive configuration
    yend = -lower_bound
    # yend = 0  # negative configuration

    # Omit the endpoints: assume 0
    y0 = y0[1: -1]

    # Sanity check
    assert np.isclose(
        length, arc(y0, x, yend), atol=1e-14, rtol=0,
    )

    result = minimize(
        fun=potential_energy, x0=y0, args=(x, yend),
        bounds=Bounds(
            lb=0, ub=yend,  # positive configuration
            # lb=lower_bound, ub=0,  # negative configuration
        ),
        constraints=NonlinearConstraint(
            fun=functools.partial(arc, x=x, yend=-lower_bound),
            lb=length, ub=length,
        ),
        options={'maxiter': 1_000},
    )
    if not result.success:
        raise ValueError(result.message)

    y, x_diff, dy_dx, integrand = unpack(x, result.x, yend)
    return x, y


def main() -> None:
    x, y = solve(resolution=49)  # 101
    plt.plot(x, y)
    plt.show()


if __name__ == '__main__':
    main()

任何高于此值都无法收敛。这可以通过提供解析雅可比行列式来帮助。