Runge-Kutta/main.ipynb

In [5]:

import numpy as np
import matplotlib.pyplot as plt

# ===================== First Task Functions ===================== #

# Equations for the first task
def func_a(x, y):
    return x + np.cos(y)

def func_b(x, y):
    return x**2 + y**2

# Runge-Kutta 2nd order method for first-order ODE
def runge_kutta_2(f, x0, y0, h, N):
    x = x0 + np.arange(N+1) * h
    y = np.zeros(N+1)
    y[0] = y0
    for k in range(N):
        k1 = f(x[k], y[k])
        k2 = f(x[k] + h, y[k] + h * k1)
        y[k+1] = y[k] + h * (k1 + k2) / 2
    return x, y

# Runge-Kutta 4th order method for first-order ODE
def runge_kutta_4(f, x0, y0, h, N):
    x = x0 + np.arange(N+1) * h
    y = np.zeros(N+1)
    y[0] = y0
    for k in range(N):
        k1 = f(x[k], y[k])
        k2 = f(x[k] + h/2, y[k] + h * k1 / 2)
        k3 = f(x[k] + h/2, y[k] + h * k2 / 2)
        k4 = f(x[k] + h, y[k] + h * k3)
        y[k+1] = y[k] + h * (k1 + 2*k2 + 2*k3 + k4) / 6
    return x, y

# ===================== Second Task Functions ===================== #

# Converts second-order ODE to a system of first-order ODEs
def second_order_to_system(x, y, dy):
    return dy, y * np.sin(x)

# Runge-Kutta 2nd order method for second-order ODE
def runge_kutta_2_second_order(f, x0, y0, dy0, h, N):
    x = x0 + np.arange(N+1) * h
    y = np.zeros(N+1)
    dy = np.zeros(N+1)
    y[0] = y0
    dy[0] = dy0
    for k in range(N):
        k1_y, k1_dy = f(x[k], y[k], dy[k])
        k2_y, k2_dy = f(x[k] + h, y[k] + h * k1_y, dy[k] + h * k1_dy)
        y[k+1] = y[k] + h * (k1_y + k2_y) / 2
        dy[k+1] = dy[k] + h * (k1_dy + k2_dy) / 2
    return x, y

# Runge-Kutta 4th order method for second-order ODE
def runge_kutta_4_second_order(f, x0, y0, dy0, h, N):
    x = x0 + np.arange(N+1) * h
    y = np.zeros(N+1)
    dy = np.zeros(N+1)
    y[0] = y0
    dy[0] = dy0
    for k in range(N):
        k1_y, k1_dy = f(x[k], y[k], dy[k])
        k2_y, k2_dy = f(x[k] + h/2, y[k] + h * k1_y / 2, dy[k] + h * k1_dy / 2)
        k3_y, k3_dy = f(x[k] + h/2, y[k] + h * k2_y / 2, dy[k] + h * k2_dy / 2)
        k4_y, k4_dy = f(x[k] + h, y[k] + h * k3_y, dy[k] + h * k3_dy)
        y[k+1] = y[k] + h * (k1_y + 2*k2_y + 2*k3_y + k4_y) / 6
        dy[k+1] = dy[k] + h * (k1_dy + 2*k2_dy + 2*k3_dy + k4_dy) / 6
    return x, y

# ===================== Main Code for Both Tasks ===================== #

# Parameters for both tasks
methods = {
    'RK2': {'first_order': runge_kutta_2, 'second_order': runge_kutta_2_second_order},
    'RK4': {'first_order': runge_kutta_4, 'second_order': runge_kutta_4_second_order}
}

h_values = [0.1, 0.05, 0.01, 0.005, 0.001]
N_values = [int((2 - 1)/h) for h in h_values]  # For the first task with x from 1 to 2

# First Task Equations
equations_first = {
    'a': {'func': func_a, 'x0': 1.0, 'y0': 30.0, 'x_end': 2.0},
    'b': {'func': func_b, 'x0': 2.0, 'y0': 1.0, 'x_end': 1.0}  # Reverse integration
}

# Second Task Parameters
x0_second = 0.0
y0_second = 0.0
dy0_second = 1.0
x_end_second = 1.0
N_values_second = [int((x_end_second - x0_second)/h) for h in h_values]

# Function to perform calculations and plotting for a given task
def solve_task(task_label, equations, N_values_task, x0_task, x_end_task, is_second_order=False):
    for eq_label, eq_params in equations.items():
        print(f"\nSolution for equation {task_label}{eq_label}:")

        x0 = eq_params['x0']
        y0 = eq_params['y0']
        if is_second_order:
            dy0 = eq_params['dy0']
        x_end = eq_params['x_end']
        if x_end < x0:
            direction = -1
        else:
            direction = 1

        for method_label, method_funcs in methods.items():
            method_func = method_funcs['second_order'] if is_second_order else method_funcs['first_order']
            print(f"\nMethod {method_label}:")

            ys = {}
            xs = {}
            abs_errors = []
            max_abs_errors = []

            for h, N in zip(h_values, N_values_task):
                h = direction * h  # Account for integration direction
                if is_second_order:
                    x, y = method_func(second_order_to_system, x0, y0, dy0, h, N)
                else:
                    x, y = method_func(eq_params['func'], x0, y0, h, N)
                xs[N] = x
                ys[N] = y

            for i in range(1, len(N_values_task)):
                N_prev = N_values_task[i-1]
                N_curr = N_values_task[i]

                # Find common indices for comparison
                factor = N_curr // N_prev
                indices_prev = [1, N_prev]
                indices_curr = [k * factor for k in indices_prev]

                # Compute relative errors
                rel_errors = np.abs(ys[N_curr][indices_curr] - ys[N_prev][indices_prev])
                print(f"N_{N_curr}: {' '.join(map(str, rel_errors))}")

                # Compute maximum absolute errors
                max_error = np.max(rel_errors)
                max_abs_errors.append(max_error)

            # Plotting log2 of absolute errors
            plt.plot(np.log2(N_values_task[1:]), np.log2(max_abs_errors), label=f"{method_label} for equation {task_label}{eq_label}")

    plt.xlabel("log2(N)")
    plt.ylabel("log2(max absolute error)")
    plt.legend()
    plt.title(f"Logarithm of absolute errors for task {task_label}")
    plt.show()

# ===================== Solve First Task ===================== #
solve_task('1', equations_first, N_values, x0_task=None, x_end_task=None)

# ===================== Solve Second Task ===================== #

# Second Task Equation
equations_second = {
    '': {'dy0': dy0_second, 'x0': x0_second, 'y0': y0_second, 'x_end': x_end_second}
}

solve_task('2', equations_second, N_values_second, x0_task=x0_second, x_end_task=x_end_second, is_second_order=True)

Solution for equation 1a:

Method RK2:
N_20: 0.0002785283564428198 0.0038035441958399474
N_100: 4.479515251887278e-05 0.0012767268215512217
N_200: 2.7750003539495083e-07 4.103739219374347e-05
N_1000: 4.440127909788316e-08 1.3191150621594261e-05

Method RK4:
N_20: 2.9636773035690567e-07 3.9579403434686355e-06
N_100: 9.184862648226044e-09 2.6419341025984977e-07
N_200: 2.6076918402395677e-12 3.9756997693984886e-10
N_1000: 8.526512829121202e-14 2.646061147970613e-11

Solution for equation 1b:

Method RK2:
N_20: 0.0014339682258658337 0.019603229003013034
N_100: 9.465097127459021e-05 0.0035429668461870456
N_200: 5.85138142383812e-08 6.395910166601126e-05
N_1000: 2.542437982366863e-08 1.8243564207320873e-05

Method RK4:
N_20: 1.6409863435429273e-05 1.5852614062783488e-05
N_100: 5.942575315165399e-07 2.292905823431113e-06
N_200: 1.8908408172535474e-10 3.99405974960132e-09
N_1000: 6.314393452555578e-12 2.6547564146994773e-10

Solution for equation 2:

Method RK2:
N_20: 6.247396158842733e-06 0.0006487724632973091
N_100: 4.999227261395789e-07 0.00021289741698327092
N_200: 6.249973958510902e-10 6.760371213720973e-06
N_1000: 4.999992173765344e-11 2.169024081455362e-06

Method RK4:
N_20: 1.778789573969597e-09 1.0195093080866968e-07
N_100: 3.0951276264179484e-11 6.3144536266435125e-09
N_200: 1.9359513991901167e-15 8.852918398360998e-12
N_1000: 3.382710778154774e-17 5.822009541134321e-13