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

def lagrangian(q, qp, m=1.0, k=1.0):
    """
    Define the Lagrangian for a simple harmonic oscillator.
        :param q: Generalized coordinate.
        :param qp: Generalized velocity.
        :param m: Mass (default is 1).
        :param k: Spring constant (default is 1).
    :return: 
        Lagrangian value.
    """
    return 0.5 * m * qp**2 - 0.5 * k * q**2


def action_discrete(q, t, m=1.0, k=1.0):
    """
    Define the Lagrangian for a simple harmonic oscillator.
        :param q: Generalized coordinate.
        :param t: Time.
        :param m: Mass (default is 1).
        :param k: Spring constant (default is 1).
    :return: 
        Action_discrete.
    """
    dt = t[1] - t[0]
    qp = np.gradient(q, dt)
    L = lagrangian(q, qp, m, k)
    return np.sum(L) * dt


def extremize_action(q_init, t1, t2, m=1.0, k=1.0):
    """
    Define the Lagrangian for a simple harmonic oscillator.
        :param q: Generalized coordinate.
        :param qp: Generalized velocity.
        :param t: Time.
        :param m: Mass (default is 1).
        :param k: Spring constant (default is 1).
    :return: 
        Lagrangian value.
    """
    return 0.5 * m * qp**2 - 0.5 * k * q**2

    return T - V 

def action_discrete(q, t, m=1.0, k=1.0):
    dt = t[1] - t[0]
    qp = np.gradient(q, dt)
    L = lagrangian(q, qp, m, k)
    return np.sum(L) * dt
    
def action_integral(q_func, t1, t2, args):
    """
    Compute the action integral for the given path q(t).
        :param q_func: Function representing the path q(t).
        :param tl: Start time.
        :param t2: End time.
        :param args: Additional arguments for the Lagrangian.
    :return: 
        Evaluated action.
    """
    integrand = lambda t: lagrangian(q_func(t), np.gradient(q_func(t)), t, *args)
    S, _ = quad(integrand, t1, t2)
    return S 

def extremize_action(q_init, t1, t2, num_points):
    """
    Numerically extremize the action using a simple gradient-based
    approach.
        :param q_init: Initial guess for the path as a function of time array.
        :param t1: Start time.
        :param t2: End time.
        :param num_points: Number of points in the discretized time domain.
    :retum: 
        Extremized path as a time array.
    """
    n = len(q_init)
    t = np.linspace(t1, t2, n)

    # boundary conditions
    q0, qN = q_init[0], q_init[-1]

    def objective(q_internal):
        q_full = np.concatenate(([q0], q_internal, [qN]))
        return action_discrete(q_full, t, m, k)

    q_internal_init = q_init[1:-1]
    result = minimize(objective, q_internal_init, method="BFGS")

    q_extremal = np.concatenate(([q0], result.x, [qN]))
    return t, q_extremal


# --- Example ---
m, k = 1, 1
initial_path = np.linspace(0.0, 1.0, 100)
t, extremal_path = extremize_action(initial_path, 0.0, 10.0, 100)
#print("Extremal Path:", extremal_path)
#print(t)
plt.plot(t, extremal_path)

[<matplotlib.lines.Line2D at 0x263b0bfef10>]

Action Functional and Its Extremals¶