Runge-Kutta/main.ipynb

In [1]:

import numpy as np
import matplotlib.pyplot as plt



# Определяем функции для уравнений
def func_a(x, y):
    return x + np.cos(y)

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

# Реализация метода Рунге-Кутты 2-го порядка
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

# Реализация метода Рунге-Кутты 4-го порядка
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

# Начальные условия и параметры
equations = {
    '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}  # Обратное интегрирование
}

methods = {
    'RK2': runge_kutta_2,
    'RK4': runge_kutta_4
}

h_values = [0.1, 0.05, 0.01, 0.005, 0.001]
N_values = [10, 20, 100, 200, 1000]

for eq_label, eq_params in equations.items():
    print(f"\nРешение для уравнения {eq_label}:")

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

    for method_label, method_func in methods.items():
        print(f"\nМетод {method_label}:")

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

        for h, N in zip(h_values, N_values):
            h = direction * h  # Учёт направления интегрирования
            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)):
            N_prev = N_values[i-1]
            N_curr = N_values[i]

            # Поиск общих индексов для сравнения
            factor = N_curr // N_prev
            indices = [int(k * factor) for k in [1, N_prev]]

            # Вычисление относительных ошибок
            rel_errors = np.abs(ys[N_curr][indices] - ys[N_prev][[1, N_prev]])
            print(f"N_{N_curr}: {' '.join(map(str, rel_errors))}")

            # Вычисление максимальных абсолютных ошибок
            max_error = np.max(rel_errors)
            max_abs_errors.append(max_error)

        # Построение графика логарифма абсолютных ошибок
        plt.plot(np.log2(N_values[1:]), np.log2(max_abs_errors), label=f"{method_label} для уравнения {eq_label}")

plt.xlabel("log2(N)")
plt.ylabel("log2(max абсолютная ошибка)")
plt.legend()
plt.title("График логарифма абсолютных ошибок")
plt.show()

Решение для уравнения a:

Метод 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

Метод 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

Решение для уравнения b:

Метод 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

Метод 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