Lab3

#include <iostream>
#include <cmath>

using std::cout;
using std::cin;
using std::endl;

static const double h = 0.1;

double derivative(double x, double y); // Дифференциальное уравнение
double correct(double x_last, double y_last, double x_next, double y_cur, double h, double eps); // Проверка предсказанного значения \\ Мод. метод Эйлера
double rungeK(double h, double x, double y, int number); // Вычисление k1, k2, k3, k4 \\ Метод Рунге-Кутта
double yt(double x, double y);
void eulerMet(bool select); // Метод Эйлера
void rungeKutt(bool select); // Метод Рунге-Кутта
void milnMet(); // Метод Милна

int main()
{
    unsigned short int pck = -1;
    bool check_task = false;
    cout << "Solution of the Cauchy problem" << endl;
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    cout << "Choose a method:\n1.Euler's method\n2.Runge-Kutt method\n3.Miln method\nEnter:";

    while (check_task == false)
    {
        cin >> pck;
        switch (pck)
        {
        case 1:
        {
            check_task = true;
            double a, b, eps;
            bool mod;
            cout << "0 - Not modified\n1 - Modified\nEnter: ";
            cin >> mod;
            if (mod)
            {
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
                cout << "Modified Euler's method" << endl;
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
            }
            else
            {
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
                cout << "Not modified Euler's method" << endl;
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
            }

            eulerMet(mod);

            break;
        }
        case 2:
        {
            check_task = true;
            double a, b, eps;
            bool mod;
            cout << "0 - Without precision\n1 - With a precision\nEnter: ";
            cin >> mod;
            if (mod)
            {
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
                cout << "Runge-Kutt method with a given precision" << endl;
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
            }
            else
            {
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
                cout << "Runge-Kutt method without precision" << endl;
                cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
            }

            rungeKutt(mod);

            break;
        }
        case 3:
        {
            cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
            cout << "\tMiln Method" << endl;
            cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;

            milnMet();

            break;
        }
        default: {cout << "#Error#\nRe-enter: "; }
        }
    }

    return 0;
}

double derivative(double x, double y)
{
    return (exp(x) + exp(-x));
}
double yt(double x)
{
    return(exp(x) - exp(-x));
}
double correct(double x_last, double y_last, double x_next, double y_cur, double h, double eps)
{
    double y_correct = y_cur;
    do
    {
        y_cur = y_correct;
        y_correct = y_last + h * 0.5 * (derivative(x_last, y_last) + derivative(x_next, y_cur));
    } while (fabs(y_correct - y_cur) > eps);

    return y_correct;
}
void eulerMet(bool select)
{
    double a = 0, b = 1;
    double n = (b - a) / h;
    double* massX = new double[n + 1];
    double* massY = new double[n + 1];
    double* massT = new double[n + 1];
    cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
    cout << "h (step) = " << h << endl;
    massX[0] = a;
    massY[0] = 0;
    switch (select)
    {
    case true:
    {
        double eps;
        double x1;
        double y_correct;
        cout << "Enter eps = "; cin >> eps;
        for (unsigned short int i = 1; i <= n; i++)
        {
            massX[i] = a + i * h;
            massY[i] = massY[i - 1] + h * derivative(massX[i - 1], massY[i - 1]);
            x1 = massX[i] + h;
            massY[i] = correct(massX[i - 1], massY[i - 1], x1, massY[i], h, eps);

        }

        break;
    }
    case false:
    {
        for (unsigned short int i = 1; i <= n; i++)
        {
            massX[i] = a + i * h;
            massY[i] = massY[i - 1] + h * derivative(massX[i - 1], massY[i - 1]);
        }

        break;
    }
    }
    for (unsigned short int i = 0; i <= n; i++)
    {
        massT[i] = yt(massX[i]);
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
    cout << "Results: " << endl;
    for (unsigned short int i = 0; i <= n; i++)
    {
        cout.precision(9);
        cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\terror = " << fabs(massY[i] - massT[i]) << endl;
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    delete[] massX;
    delete[] massY;
}
void rungeKutt(bool select)
{
    double a = 0.0, b = 1.0;
    double n = (b - a) / h;
    double* massX = new double[n + 1];
    double* massY = new double[n + 1];
    double* massT = new double[n + 1];
    double k1, k2, k3, k4;
    cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
    cout << "h (step) = " << h << endl;
    massX[0] = a;
    massY[0] = 0;
    switch (select)
    {
    case true:
    {
        double eps;
        double h0 = (b - a) / n;
        double h1 = h0 / 2.0;
        double y1, y2, yv;
        double sum1 = 0, sum2 = 0, sum3 = 0;
        cout << "Enter eps = "; cin >> eps;
        for (unsigned short int i = 1; i <= n; i++)
        {
            massX[i] = a + i * h0;
            for (unsigned short int j = 1; i < 5; i++)
            {
                sum1 += h0 * rungeK(h0, massX[i - 1], massY[i - 1], j);
            }
            y1 = massY[i - 1] + sum1;
            for (unsigned short int j = 1; i < 5; i++)
            {
                sum2 += h1 * rungeK(h1, massX[i - 1], massY[i - 1], j);
            }
            yv = massY[i - 1] + sum2;
            for (unsigned short int j = 1; i < 5; i++)
            {
                sum3 += h1 * rungeK(h1, massX[i - 1], yv, j);
            }
            y2 = yv + sum3;

            if (fabs(y1 - y2) < eps)
            {
                for (unsigned short int i = 1; i <= n; i++)
                {
                    massX[i] = a + i * h0;
                    massY[i] = massY[i - 1] + (((rungeK(h0, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h0, massX[i - 1],
                        massY[i - 1], 2) + rungeK(h0, massX[i - 1], massY[i - 1], 3)) + rungeK(h0, massX[i - 1], massY[i - 1], 4)) / 6.0);
                }
            }
            else
            {
                for (unsigned short int i = 1; i <= n; i++)
                {
                    massX[i] = a + i * h0;
                    massY[i] = massY[i - 1] + (((rungeK(h1, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h1, massX[i - 1],
                        massY[i - 1], 2) + rungeK(h1, massX[i - 1], massY[i - 1], 3)) + rungeK(h1, massX[i - 1], massY[i - 1], 4)) / 6.0);
                }
            }
        }
        break;
    }
    case false:
    {
        for (unsigned short int i = 1; i <= n; i++)
        {
            massX[i] = a + i * h;
            massY[i] = massY[i - 1] + (((rungeK(h, massX[i - 1], massY[i - 1], 1)) + 2 * (rungeK(h, massX[i - 1],
                massY[i - 1], 2) + rungeK(h, massX[i - 1], massY[i - 1], 3)) + rungeK(h, massX[i - 1], massY[i - 1], 4)) / 6.0);
        }
        break;
    }
    }
    for (unsigned short int i = 0; i <= n; i++)
    {
        massT[i] = yt(massX[i]);
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
    cout << "Results: " << endl;
    for (unsigned short int i = 0; i <= n; i++)
    {
        cout.precision(9);
        cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\t"; printf("error = % 2.9f\n", fabs(massY[i] - massT[i]));
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    delete[] massX;
    delete[] massY;
}
void milnMet()
{
    double a = 0.0, b = 1.0;
    double n = (b - a) / h;
    double* massX = new double[n + 1];
    double* massY = new double[n + 1];
    double* massT = new double[n + 1];
    cout << "a (left border) = " << a << "\nb (right border) = " << b << endl;
    cout << "h (step) = " << h << endl;
    massX[0] = a;
    for (unsigned short int i = 1; i <= n; i++)
        massX[i] = a + i * h;

    massY[0] = 0.0; massY[1] = 0.200333507; massY[2] = 0.402672019; massY[3] = 0.609040608;
    double temp;
    for (unsigned short int i = 4; i <= n; i++)
    {
        temp = massY[i - 4] + (4 * h * (2 * derivative(massX[i - 3], massY[i - 3]) - derivative(massX[i - 2], massY[i - 2])
            + 2 * derivative(massX[i - 1], massY[i - 1])) / 3.0);
        massY[i] = massY[i - 2] + (h * (derivative(massX[i - 2], massY[i - 2]) + 4 * derivative(massX[i - 1], massY[i - 1]) +
            derivative(massX[i], temp)) / 3.0);
    }
    for (unsigned short int i = 0; i <= n; i++)
    {
        massT[i] = yt(massX[i]);
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;
    cout << "Equation: y` = e^x + e^(-x); y(0) = 0; yt = e^x - e^(-x)" << endl;
    cout << "Results: " << endl;
    for (unsigned short int i = 0; i <= n; i++)
    {
        cout.precision(9);
        cout << "X[" << i << "] = " << massX[i] << "\tY[" << i << "] = " << massY[i] << "\t"; printf("error = % 2.9f\n", fabs(massY[i] - massT[i]));
    }
    cout << "~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~" << endl;

    delete[] massX;
    delete[] massY;
    delete[] massT;
}
double rungeK(double h, double x, double y, int number)
{
    double result1, result2, result3, result4;
    switch (number)
    {
    case 1:
    {
        result1 = h * derivative(x, y);
        return result1;
        break;
    }
    case 2:
    {
        result1 = h * derivative(x, y);
        result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
        return result2;
        break;
    }
    case 3:
    {
        result1 = h * derivative(x, y);
        result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
        result3 = h * derivative(x + (h / 2.0), y + (result2 / 2.0));
        return result3;
        break;
    }
    case 4:
    {
        result1 = h * derivative(x, y);
        result2 = h * derivative(x + (h / 2.0), y + (result1 / 2.0));
        result3 = h * derivative(x + (h / 2.0), y + (result2 / 2.0));
        result4 = h * derivative(x + h, y + result3);
        return result4;
        break;
    }
    default:cout << "ERROR_K" << endl;
    }
}