added Runge Kutta

runge-kutta.py

@@ -0,0 +1,78 @@
+import numpy as np
+
+def runge_kutta(n, y, f, x, h, m, k):
+    """
+    Perform Runge-Kutta calculation using Gill's method.
+
+    Args:
+        n (int): Number of differential equations to be solved.
+        y (np.ndarray): Array of N dependent variables with initial values.
+        f (np.ndarray): Array of the N derivatives of the variables Y.
+        x (float): Independent variable.
+        h (float): Step size.
+        m (int): Index used in the subroutine.
+        k (int): Integer used as the argument of a computed GO TO statement.
+
+    Returns:
+        Tuple[int, int]: Updated values of m and k.
+    """
+    m += 1
+
+    if m == 1:
+        q = np.zeros(n)
+        a = 0.5
+    elif m == 2:
+        a = 1.707107
+        x += 0.5 * h
+        q = np.zeros(n)
+        for i in range(n):
+            y[i] = y[i] + a * (f[i] * h - q[i])
+            q[i] = 2.0 * a * h * f[i] + (1.0 - 3.0 * a) * q[i]
+        a = 0.2928932
+    elif m == 3:
+        a = 1.707107
+    elif m == 4:
+        x += 5 * h
+        a = 0.2928932
+
+    if a == 0.2928932:
+        for i in range(n):
+            y[i] = y[i] + h * f[i] / 6.0 - q[i] / 3.0
+        m = 0
+        k = 2
+
+    return m, k
+
+def main():
+    """
+    Main program to solve differential equations using Runge-Kutta method.
+
+    User provides the initial values for the differential equation, and the
+    program calculates the solutions.
+    """
+    n =  # Set the number of differential equations
+    y = np.zeros(n)
+    f = np.zeros(n)
+    q = np.zeros(n)
+
+    x, n_eq, xlim, h, m, *y_values = map(float, input("Enter initial values: ").split())
+    y[:n_eq] = y_values
+
+    k = 1
+    while True:
+        if x > xlim:
+            break
+
+        if k == 1:
+            # Programmer must provide expressions for each F() in the main program
+            f[0] =  # Add the expression for the 1st derivative
+            # Add similar lines for other derivatives if n > 1
+            print(f"At X={x}, Y={y[0]}")
+        elif k == 2:
+            k = 1
+            m, k = runge_kutta(n_eq, y, f, x, h, m, k)
+
+        x += h
+
+if __name__ == "__main__":
+    main()