Wronskian for second order differential equations

maths/wronskian_second_order_de.py

@@ -0,0 +1,143 @@
+"""
+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
+