Lab 6: Solution of Systems of ODEs using Runge-Kutta 4th Order Method

import numpy as np import matplotlib.pyplot as plt def F(x, y, z): # y’ return 3*y+4*z def G(x, y, z): # z’ return 3*z-4*y def phi1(x): return np.exp(3*x)*np.cos(4*x) def phi2(x): return -np.exp(3*x)*np.sin(4*x) def RK4(f, g, X, y0, z0): n = len(X) Y = np.zeros(n, dtype=float) Z = np.zeros(n, dtype=float) Y[0] = y0 Z[0] = z0 for i in range(n – 1): h = X[i + 1] – X[i] # First set of weights k1 = h * f(X[i], Y[i], Z[i]) l1 = h * g(X[i], Y[i], Z[i]) # Second set of weights k2 = h * f(X[i] + h/2, Y[i] + k1/2, Z[i] + l1/2) l2 = h * g(X[i] + h/2, Y[i] + k1/2, Z[i] + l1/2) # Third set of weights k3 = h * f(X[i] + h/2, Y[i] + k2/2, Z[i] + l2/2) l3 = h * g(X[i] + h/2, Y[i] + k2/2, Z[i] + l2/2) # Fourth set of weights k4 = h * f(X[i] + h, Y[i] + k3, Z[i] + l3) l4 = h * g(X[i] + h, Y[i] + k3, Z[i] + l3) # Weighted average k = (k1 + 2*k2 + 2*k3 + k4) / 6 l = (l1 + 2*l2 + 2*l3 + l4) / 6 Y[i + 1] = Y[i] + k Z[i + 1] = Z[i] + l return Y, Z # Initial conditions x0, y0, z0 = 0.0, 1.0, 0.0 xn, n = 1, 20 X = np.linspace(x0, xn, n+1) # Solve using RK4 yRK4, zRK4 = RK4(F, G, X, y0, z0) # Exact solutions yExact = phi1(X) zExact = phi2(X) # Print results print(f”{‘i’: < 3}{'x': < 10}{'yRK4': < 10}{'yExact':< 10}{'zRK4': < 10}{'zExact':< 10}") print("-" * 65) for i in range(n + 1): print(f"{i:< 3}{X[i]: < 10.5f}{yRK4[i]: < 10.5f}{yExact[i]: < 10.5f}{zRK4[i]: < 10.5f}{zExact[i]: < 10.5f}") # Plot results plt.figure(figsize=(10, 6)) plt.plot(X, yExact, label="Exact y(x)", color="black", marker="o") plt.plot(X, yRK4, label="RK4 y(x)", color="m", marker="d") plt.plot(X, zExact, label="Exact z(x)", color="b", linestyle='--', marker="o") plt.plot(X, zRK4, label="RK4 z(x)", color="m", linestyle='--', marker="d") plt.title("Comparison of Numerical Methods for System of ODEs") plt.xlabel("x", fontsize=12) plt.ylabel("y, z", fontsize=12) plt.legend(fontsize=10, edgecolor="black") plt.grid(True, linestyle="--") plt.show()