wronskian_second_order_de.py

"""
wronskian_second_order_de.py
A symbolic and numerical exploration of the Wronskian for second-order linear differential equations.

This program:
1. Takes coefficients (a, b, c) for a*y'' + b*y' + c*y = 0.
2. Computes characteristic roots.
3. Classifies the solution type.
4. Constructs the general solution.
5. Demonstrates Wronskian computation.

Author: Venkat Thadi

References: https://tutorial.math.lamar.edu/classes/de/wronskian.aspx
"""

import math
import cmath

def compute_characteristic_roots(a: float, b: float, c: float) -> tuple[complex, complex]:
    """
    Compute characteristic roots for a second-order homogeneous linear DE.

    >>> compute_characteristic_roots(1, -3, 2)
    (2.0, 1.0)
    >>> compute_characteristic_roots(1, 2, 5)
    ((-1+2j), (-1-2j))
    """
    if a == 0:
        raise ValueError("Coefficient 'a' cannot be zero for a second-order equation.")

    discriminant = b ** 2 - 4 * a * c
    sqrt_disc = cmath.sqrt(discriminant)
    root1 = (-b + sqrt_disc) / (2 * a)
    root2 = (-b - sqrt_disc) / (2 * a)

    # Simplify if roots are purely real
    if abs(root1.imag) < 1e-12:
        root1 = float(root1.real)
    if abs(root2.imag) < 1e-12:
        root2 = float(root2.real)

    return root1, root2


def classify_solution_type(root1: complex, root2: complex) -> str:
    """
    Classify the nature of the roots.

    >>> classify_solution_type(2, 1)
    'Distinct Real Roots'
    >>> classify_solution_type(1+2j, 1-2j)
    'Complex Conjugate Roots'
    >>> classify_solution_type(3, 3)
    'Repeated Real Roots'
    """
    if isinstance(root1, complex) and isinstance(root2, complex) and root1.imag != 0:
        return "Complex Conjugate Roots"
    elif root1 == root2:
        return "Repeated Real Roots"
    else:
        return "Distinct Real Roots"


def compute_wronskian(f, g, f_prime, g_prime, x: float) -> float:
    """
    Compute Wronskian determinant W(f, g) = f * g' - f' * g.

    >>> import math
    >>> def f(x): return math.exp(x)
    >>> def g(x): return x * math.exp(x)
    >>> def f_prime(x): return math.exp(x)
    >>> def g_prime(x): return math.exp(x) + x * math.exp(x)
    >>> round(compute_wronskian(f, g, f_prime, g_prime, 0), 3)
    1.0
    """
    return f(x) * g_prime(x) - f_prime(x) * g(x)


def construct_general_solution(root1: complex, root2: complex) -> str:
    """
    Construct the general solution based on the roots.

    >>> construct_general_solution(2, 1)
    'y(x) = C1 * e^(2x) + C2 * e^(1x)'
    >>> construct_general_solution(3, 3)
    'y(x) = (C1 + C2 * x) * e^(3x)'
    >>> construct_general_solution(-1+2j, -1-2j)
    'y(x) = e^(-1x) * (C1 * cos(2x) + C2 * sin(2x))'
    """
    if isinstance(root1, complex) and root1.imag != 0:
        alpha = round(root1.real, 10)
        beta = round(abs(root1.imag), 10)
        return f"y(x) = e^({alpha:g}x) * (C1 * cos({beta:g}x) + C2 * sin({beta:g}x))"
    elif root1 == root2:
        return f"y(x) = (C1 + C2 * x) * e^({root1:g}x)"
    else:
        return f"y(x) = C1 * e^({root1:g}x) + C2 * e^({root2:g}x)"


def analyze_differential_equation(a: float, b: float, c: float) -> None:
    """
    Analyze the DE and print the roots, type, and general solution.

    >>> analyze_differential_equation(1, -3, 2)  # doctest: +ELLIPSIS
    Characteristic Roots: (2.0, 1.0)
    Solution Type: Distinct Real Roots
    General Solution: y(x) = C1 * e^(2x) + C2 * e^(1x)
    """
    roots = compute_characteristic_roots(a, b, c)
    root1, root2 = roots
    sol_type = classify_solution_type(root1, root2)
    general_solution = construct_general_solution(root1, root2)

    print(f"Characteristic Roots: ({root1:.1f}, {root2:.1f})")
    print(f"Solution Type: {sol_type}")
    print(f"General Solution: {general_solution}")


def main() -> None:
    """
    Main function to run the second-order differential equation Wronskian analysis.

    Interactive input is expected, so this function is skipped in doctests.
    """
    print("Enter coefficients for the equation a*y'' + b*y' + c*y = 0")

    # Skipping main in doctests because input() cannot be tested directly
    a = float(input("a = ").strip())  # doctest: +SKIP
    b = float(input("b = ").strip())  # doctest: +SKIP
    c = float(input("c = ").strip())  # doctest: +SKIP

    if a == 0:
        print("Invalid input: coefficient 'a' cannot be zero.")
        return

    analyze_differential_equation(a, b, c)



if __name__ == "__main__":
    main() # doctest: +SKIP