RK4

lab4.cpp

#include <iostream>
#include <cstdint>
#include <cmath>
#include <vector>

/*
 * y''(x) + p(x) * y'(x) + q(x) * y(x) = f(x)
 * y''(x) + y'(x) - y(x) = f(x)
 * y(x) = e^x, p(x) = 1, q(x) = -1
 * Начальные условия з. Коши:
 * у(0) = 1, y'(0) = 1, x \in [0, 1]
*/


const float a = 0.0f, b = 1.0f;
const uint32_t n = 10;
const float h = (b - a) / n;

inline static float
RK_iterate(float y_n, const std::vector<float> &k1234)
{
    return(y_n + h / 6.0f * (k1234[0] + 2.0f * k1234[1] + 2.0f * k1234[2] + k1234[3]));
}

static std::vector<float>
RK_1234(float x_n)
{
    float k_1 = std::exp(x_n);
    float k_2 = std::exp(x_n + h / 2.0f) - h / 2.0f * k_1;
    float k_3 = std::exp(x_n + h / 2.0f) - h / 2.0f * k_2;
    float k_4 = std::exp(x_n + h) - h * k_3;

    return {k_1, k_2, k_3, k_4};
}

static float
RK(float y_0)
{
    float y_n = y_0;
    for (uint32_t i = 0; i <= n; ++i)
    {
        auto k1234 = RK_1234(h * n);
        y_n = RK_iterate(y_n, k1234);
    }

    return(y_n);
}

int32_t
main(void)
{
    std::cout << RK(1.0f) << std::endl;
    return(0);
}