Runge-Kutta/main.ipynb

In [29]:

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]
legend_abbreviation = "EMSH64D04"

# 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]

                indices_prev = np.linspace(0, N_prev, 11, endpoint=True, dtype=int)
                indices_curr = np.linspace(0, N_curr, 11, endpoint=True, dtype=int)
                
                rel_errors = np.abs(ys[N_curr][indices_curr] - ys[N_prev][indices_prev])

                rel_errors_without_zero = rel_errors[1:]

                if len(rel_errors_without_zero) != 10:
                    raise ValueError(f"The resulting error array does not contain exactly 10 values (found {len(rel_errors_without_zero)}).")

                print(f"N_{N_curr}:")
                print(" ".join(f"{val:.10g}" for val in rel_errors_without_zero))

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

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


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

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

# ===================== 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.0002785283564 0.0006147262836 0.001008931808 0.001453439014 0.001930378282 0.002412818397 0.002869546941 0.003270516965 0.003589271118 0.003803544196
N_100:
9.442991941e-05 0.0002090095219 0.0003440780151 0.0004970853173 0.000661623814 0.0008275460164 0.000982606614 0.001114948045 0.001215060742 0.001276726822
N_200:
3.056533167e-06 6.776709096e-06 1.117548006e-05 1.617051613e-05 2.154595543e-05 2.695219014e-05 3.196167493e-05 3.616337836e-05 3.92460027e-05 4.103739219e-05
N_1000:
9.836686345e-07 2.181501333e-06 3.598507423e-06 5.208179537e-06 6.94057232e-06 8.682048339e-06 1.029342858e-05 1.164114701e-05 1.262499773e-05 1.319115062e-05

Method RK4:
N_20:
2.963677304e-07 7.195955476e-07 1.286779597e-06 1.981418709e-06 2.730011495e-06 3.402355173e-06 3.863917016e-06 4.06174507e-06 4.057771658e-06 3.957940343e-06
N_100:
2.028602708e-08 4.911409945e-08 8.745490021e-08 1.338870206e-07 1.831774981e-07 2.266869501e-07 2.56161016e-07 2.690227525e-07 2.696008536e-07 2.641934103e-07
N_200:
3.11146664e-11 7.517186873e-11 1.334505839e-10 2.034639124e-10 2.769553475e-10 3.409823535e-10 3.83870713e-10 4.027747025e-10 4.044338198e-10 3.975699769e-10
N_1000:
2.071232075e-12 5.012879001e-12 8.885336911e-12 1.354294454e-11 1.841016228e-11 2.266986598e-11 2.55298005e-11 2.677325028e-11 2.689759526e-11 2.646061148e-11

Solution for equation 1b:

Method RK2:
N_20:
0.001433968226 0.002966446899 0.004500158049 0.006118077901 0.007943671725 0.01011452626 0.01276400597 0.01593449952 0.01917439382 0.019603229
N_100:
0.0002544574831 0.0006490274359 0.001083585901 0.001555473082 0.002087532901 0.002708169175 0.003435162915 0.004225851337 0.004773274812 0.003542966846
N_200:
4.989075659e-06 1.577579261e-05 2.841768018e-05 4.243009511e-05 5.827699515e-05 7.659969145e-05 9.754117317e-05 0.0001187768894 0.0001276600829 6.395910167e-05
N_1000:
1.465472705e-06 4.844589588e-06 8.845072874e-06 1.329479324e-05 1.833103302e-05 2.414770351e-05 3.077111962e-05 3.740934404e-05 3.986381442e-05 1.824356421e-05

Method RK4:
N_20:
1.640986344e-05 1.927399391e-05 1.896334616e-05 1.858634685e-05 1.928175228e-05 2.168188937e-05 2.624829711e-05 3.257189624e-05 3.426085818e-05 1.585261406e-05
N_100:
8.782501411e-07 1.029525818e-06 9.941156065e-07 9.473647006e-07 9.576468701e-07 1.064105619e-06 1.294298038e-06 1.614401477e-06 1.566632217e-06 2.292905823e-06
N_200:
1.090825874e-09 1.270972577e-09 1.199204319e-09 1.104448866e-09 1.081153056e-09 1.186212129e-09 1.462583832e-09 1.869790101e-09 1.756605084e-09 3.99405975e-09
N_1000:
7.065736884e-11 8.222408865e-11 7.727574136e-11 7.075251496e-11 6.887024284e-11 7.540712499e-11 9.326694972e-11 1.199462751e-10 1.128612759e-10 2.654756415e-10

Solution for equation 2:

Method RK2:
N_20:
6.247396159e-06 2.487878023e-05 5.565768745e-05 9.842805737e-05 0.0001533006145 0.000220811222 0.0003020454697 0.0003987265304 0.0005132652519 0.0006487724633
N_100:
1.997366387e-06 7.95941899e-06 1.783245443e-05 3.160261921e-05 4.935410358e-05 7.131823078e-05 9.791179558e-05 0.0001297638108 0.0001677303501 0.000212897417
N_200:
6.24108064e-08 2.488453061e-07 5.580570304e-07 9.90315554e-07 1.549198012e-06 2.243093594e-06 3.08637305e-06 4.100198781e-06 5.312965976e-06 6.760371214e-06
N_1000:
1.99718817e-08 7.963985801e-08 1.786277719e-07 3.170586584e-07 4.961266403e-07 7.185768469e-07 9.890781877e-07 1.314481255e-06 1.703958302e-06 2.169024081e-06

Method RK4:
N_20:
1.778789574e-09 4.842148099e-09 7.349121756e-09 7.663824464e-09 4.391501718e-09 3.636483226e-09 1.741321887e-08 3.782447944e-08 6.571693723e-08 1.019509308e-07
N_100:
1.05174397e-10 3.19153648e-10 5.187487662e-10 5.928556535e-10 4.445525059e-10 1.001532191e-11 8.449324573e-10 2.1288179e-09 3.928984693e-09 6.314453627e-09
N_200:
1.59774971e-13 5.087319455e-13 8.605893775e-13 1.045441511e-12 9.127143485e-13 3.29625216e-13 8.243405958e-13 2.663091969e-12 5.300426764e-12 8.852918398e-12
N_1000:
1.076916334e-14 3.433364704e-14 5.845324225e-14 7.138734048e-14 6.32827124e-14 2.553512957e-14 5.051514762e-14 1.717515019e-13 3.466116283e-13 5.822009541e-13