#!/usr/bin/env python3.10 """Approximate a buckling-load root with a safeguarded Newton method.""" from pathlib import Path import shutil import matplotlib as mpl import matplotlib.pyplot as plt import numpy as np def evaluate_load_function(area, stress, eccentricity, factor, radius, length, modulus, load_value): """ Evaluate the nonlinear load equation. Parameters ---------- area, stress, eccentricity, factor, radius, length, modulus : float Model parameters. load_value : float or numpy.ndarray Load value or array of load values. """ load_value = np.asarray(load_value, dtype=float) complex_argument = np.sqrt(load_value.astype(complex) / (area * modulus)) secant_term = 1.0 / np.cos((length / (2.0 * radius)) * complex_argument) function_values = load_value / area + stress / ( 1.0 + (eccentricity * factor) / (radius**2 * secant_term) ) return np.real_if_close(function_values, tol=1000).real def evaluate_load_derivative(area, stress, eccentricity, factor, radius, length, modulus, load_value): """ Approximate the derivative of the nonlinear load equation. Parameters ---------- area, stress, eccentricity, factor, radius, length, modulus : float Model parameters. load_value : float Evaluation point. """ step_size = max(1.0e-6, 1.0e-6 * abs(load_value)) return ( evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, load_value + step_size ) - evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, load_value - step_size ) ) / (2.0 * step_size) def find_sign_change_interval(area, stress, eccentricity, factor, radius, length, modulus): """ Locate a sign-change interval for the real root of interest. Returns ------- tuple[float, float] Interval where the function changes sign. """ grid_points = np.linspace(-50000.0, 0.0, 4000) function_values = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, grid_points ) for point_index in range(1, len(grid_points)): if function_values[point_index - 1] * function_values[point_index] < 0.0: return grid_points[point_index - 1], grid_points[point_index] raise RuntimeError("Unable to bracket the root.") def compute_critical_load(area, modulus, length, radius): """Compute the Euler critical load.""" return (np.pi**2 * modulus * area) / ((length / radius) ** 2) def solve_newton_method( area, stress, eccentricity, factor, radius, length, modulus, bracket_interval, initial_guess, tolerance, max_iterations, ): """ Solve the nonlinear load equation with a safeguarded Newton method. Returns ------- tuple[float, int] Approximated root and iteration count. """ left_bound, right_bound = bracket_interval current_load = min(max(float(initial_guess), left_bound), right_bound) left_value = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, left_bound ) for iteration in range(1, max_iterations + 1): derivative_value = evaluate_load_derivative( area, stress, eccentricity, factor, radius, length, modulus, current_load ) function_value = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, current_load ) if derivative_value == 0.0 or not np.isfinite(derivative_value): next_load = 0.5 * (left_bound + right_bound) else: next_load = current_load - function_value / derivative_value if next_load < left_bound or next_load > right_bound or not np.isfinite(next_load): next_load = 0.5 * (left_bound + right_bound) next_value = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, next_load ) if left_value * next_value <= 0.0: right_bound = next_load else: left_bound = next_load left_value = next_value if abs(next_load - current_load) <= tolerance or abs(next_value) <= tolerance: return next_load, iteration current_load = next_load return current_load, max_iterations def configure_plot_style(): """Configure a consistent LaTeX-like plot style.""" mpl.rcParams.update( { "text.usetex": shutil.which("latex") is not None, "font.family": "serif", "font.serif": ["Computer Modern Roman", "DejaVu Serif"], "mathtext.fontset": "cm", "axes.titlesize": 15, "axes.labelsize": 13, "xtick.labelsize": 11, "ytick.labelsize": 11, "legend.fontsize": 11, } ) def main(): """Run the buckling-load experiment and save the resulting figure.""" configure_plot_style() area = 1.0 stress = 40000.0 factor = 1.0 radius = 1.0 modulus = 30000000.0 eccentricity = 0.3 length = 50.0 critical_load = compute_critical_load(area, modulus, length, radius) bracket_interval = find_sign_change_interval( area, stress, eccentricity, factor, radius, length, modulus ) root_value, iteration_count = solve_newton_method( area, stress, eccentricity, factor, radius, length, modulus, bracket_interval, initial_guess=critical_load / 2.0, tolerance=1.0e-5, max_iterations=200, ) root_function_value = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, root_value ) load_grid = np.linspace(-50000.0, critical_load + 100.0, 1500) function_values = evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, load_grid ) output_directory = Path(__file__).resolve().parent.parent / "Figuras" output_directory.mkdir(exist_ok=True) plt.figure(figsize=(7.0, 5.0)) plt.plot(load_grid, function_values, "-k", linewidth=2.0, label="Load function") plt.axhline(0.0, color="tab:blue", linewidth=1.5, label="Zero reference") plt.plot(root_value, root_function_value, "or", label="Newton root") plt.plot( critical_load, evaluate_load_function( area, stress, eccentricity, factor, radius, length, modulus, critical_load ), "x", color="tab:red", markersize=8, label=r"$P_{\mathrm{cr}}$", ) plt.xlabel("Load") plt.ylabel("Function value") plt.legend() plt.grid(True, linewidth=0.4, alpha=0.5) plt.tight_layout() plt.savefig(output_directory / "buckling_load.pdf", bbox_inches="tight") plt.close() print(f"critical_load = {critical_load:.10f}") print(f"sign_change_interval = ({bracket_interval[0]:.10f}, {bracket_interval[1]:.10f})") print(f"root_value = {root_value:.10f}") print(f"root_function_value = {root_function_value:.10e}") print(f"iterations = {iteration_count}") if __name__ == "__main__": main()