NumericalMethods

//NumericalMethods.h

#ifndef NUMERICAL_METHODS
#define NUMERICAL_METHODS

void NotEnoughMemory();
inline double function(double x, double y);
double ** CreateMassA(const int SLAR, double **massA);
double * CreateMassB(const int SLAR, double *massB);
void CreateMass1(const int SLAR, double **mass);
double ** CreateMass2(const int SLAR, double **mass);
void swap (double *a, double *b);
void temp (const int SLAR, double **massA, double **tempA, const double *massB, double * tempB);
void PrintMass(const int SLAR, double **massA, const double *massB);
void PrintLU (const int SLAR, double **massL, double **massU);
void PrintX(const int SLAR, const double *massX);
void Gauss(const int SLAR, double **massA, double *massB);
void LU(const int SLAR, double **massA, double *massB);
void Jacobi(const int SLAR, double **massA, double *massB, double EPS);
void Zeidel(const int SLAR, double **massA, double *massB, double EPS);
double InterpolateLagrangePolynomial (const double x, double* X, double* Y, int size);
double InterpolateNewtonPolynomial(const double x, double* X, double* Y, const int size);
void MinSquare(int iPoints, double* X, double* Y, int iPolinom);
void Euler(const double x0, const double y0, const double x1, const double h);
void Runge_Kutta(const double x0, const double y0, const double x1, const double h);
void Adams(const double x0, const double y0, const double x, const double h);
#endif


//NumericalMethods.cpp

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <math.h>
#include <conio.h>
#include "NumericalMethods.h"

using namespace std;

void NotEnoughMemory()
{
    printf ("ERROR!!! Not enough memory!\n");
    system ("pause");
    exit (1);
}

inline double function(double x, double y)
{
    return x + cos(y/sqrt(5));
}

double ** CreateMassA(const int SLAR, double **massA)
{
    massA = (double **) malloc (SLAR*sizeof (double *));
    if (!massA)
    {
        NotEnoughMemory();
    }
    for (int i = 0; i < SLAR; ++i)
    {
        massA[i] = (double *)malloc (SLAR*sizeof (double));
        if (!massA[i])
        {
            NotEnoughMemory();
        }
        for (int j = 0; j < SLAR; ++j)
        {
            printf ("Enter the a[%i][%i]: ", i+1, j+1);
            scanf_s ("%lf", &massA[i][j]);
        }
    }
    return massA;
}

double * CreateMassB(const int SLAR, double *massB)
{
    massB = (double *) malloc (SLAR*sizeof (double));
    if(!massB)
    {
        NotEnoughMemory();
    }
    for (int i = 0; i < SLAR; ++i)
    {
        printf ("Enter the b[%d]: ", i + 1);
        scanf_s ("%lf", &massB[i]);
    }
    return massB;
}

double ** CreateMass2(const int SLAR, double **mass)
{
    mass = (double **) malloc (SLAR*sizeof (double *));
    if (!mass)
    {
        NotEnoughMemory();
    }
    for (int i = 0; i < SLAR; ++i)
    {
        mass[i] = (double *)malloc (SLAR*sizeof (double));
        if (!mass[i])
        {
            NotEnoughMemory();
        }
    }
    return mass;
}

void CreateMass1(const int SLAR, double **mass)
{
    *mass = (double *) malloc (SLAR*sizeof (double));
    if(!mass)
    {
        NotEnoughMemory();
    }
}

void swap (double *a, double *b)
{
    double temp = *a;
    *a = *b;
    *b = temp;
}

void temp (const int SLAR, double **massA, double **tempA, const double *massB, double * tempB)
{
    for (int i = 0; i < SLAR; ++i)
    {
        for (int j = 0; j < SLAR; ++j)
            tempA[i][j] = massA[i][j];
        tempB[i] = massB[i];
    }
}

void PrintMass(const int SLAR, double **massA, const double *massB)
{
    int k = (SLAR+1)*8+5;
    for (int i = 0; i < k; ++i)
        putchar ('=');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        putchar ('(');
        for (int j = 0; j < SLAR; ++j)
            printf (" %7.2lf", massA[i][j]);
        printf (" | %7.2lf )\n", massB[i]);
    }
}

void PrintLU (const int SLAR, double **massL, double **massU)
{
    int k = SLAR*8+6;
    int m = SLAR/2;
    for (int i = 0; i < k; ++i)
        putchar ('=');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        if (i == m) printf ("L = (");
        else printf ("    (");
        for (int j = 0; j < SLAR; ++j)
            printf (" %7.2lf", massL[i][j]);
        putchar(')'); putchar('\n');
    }
    for (int i = 0; i < k; ++i)
        putchar ('-');
    putchar('\n');
    for (int i = 0; i < SLAR; ++i)
    {
        if (i == m) printf ("U = (");
        else printf ("    (");
        for (int j = 0; j < SLAR; ++j)
            printf (" %7.2lf", massU[i][j]);
        putchar(')'); putchar('\n');
    }
}

void PrintX(const int SLAR, const double *massX)
{
    putchar ('\n');
    for (int i = 0; i < SLAR; ++i)
        printf ("x[%d] = %+.10lf\n", i+1, massX[i]);
    putchar ('\n');
}

void Gauss(const int SLAR, double **massA, double *massB)
{
    double *massX;
    CreateMass1(SLAR, &massX);
    printf ("Press 1 if you want to see step output of method of Gauss\n");
    char ch = _getch();
    if (ch == '1') PrintMass(SLAR, massA, massB);
    for (int k = 0; k < SLAR - 1; ++k)  //k - номер рівняння, над яким працюємо
    {
        int check = 0;
        double MAX = fabs(massA[k][k]);
        for (int i = k + 1; i < SLAR; ++i)  //і - номер рівняння, що іде після k-го рівняння
            if (MAX < fabs(massA[i][k]))    //визначаємо найбільше за модулем число у стовпці
            {
                check = 1;
                for (int a = 0; a < SLAR; ++a) //а - номер елемента в рівнянні
                    swap (&massA[k][a], &massA[i][a]);
                swap (&massB[k], &massB[i]);
            }
        if (MAX == 0)
        {
            for (int j = k; j < SLAR-1; ++j)
            {
                int check = 0;
                for (int i = 0; i < SLAR - 1; ++i)
                {
                    if (massA[j][i] != 0)
                    {
                        check = 1; continue;
                    }
                }
                if (massB[j] != 0)
                {
                    printf ("\nERROR!  System is indefinite\n", k+1);
                    system ("pause");
                    exit(1);
                }
            }
            printf ("\nERROR!  System is compatible\n", k+1);
            system ("pause");
            exit(1);
        }
        if ((ch == '1') && (check == 1)) PrintMass(SLAR, massA, massB);
        for (int j = k + 1; j < SLAR; ++j) //j - номер рівняння, над яким працюємо
        {                                  //занулюємо елементи під головною діагоналлю
            double a = massA[j][k] / massA[k][k]; // a - число, на яке ділиться рівняння
            for (int n = k; n < SLAR; ++n) //n - номер елемента в рівнянні
            {
                massA[j][n] = massA[j][n]-a*massA[k][n];
            }
            massB[j] = massB[j] - a*massB[k];
        }
        if (ch == '1') PrintMass(SLAR, massA, massB);
    }

    if (massA[SLAR-1][SLAR-1] == 0)
    {
        if (massB[SLAR-1] != 0) printf ("\nERROR! System is indefinite!\n");
        else printf ("\nERROR! System is compatible!\n");
        system ("pause");
        exit (1);
    }

    for (int i = SLAR-1; i >= 0; --i) //визначаємо корені рівняння
    {
        double x = 0;
        for (int j = SLAR-1; j > i; --j)
            x = x + massA[i][j]*massX[j];
        massX[i] = (massB[i] - x) / massA[i][i];
    }
    PrintX(SLAR, massX);
    free (massX);
    massX = NULL;
}

void LU(const int SLAR, double **massA, double *massB)
{
    double *massX;
    double **massL = CreateMass2(SLAR, massL);
    double **massU = CreateMass2(SLAR, massU);
    CreateMass1(SLAR, &massX);
    for (int i = 0; i < SLAR; ++i)
        for (int j = 0; j < SLAR; ++j)
        {
            if (i == j) {massU[i][j] = 1; massL[i][j] = 0;}
            else {massL[i][j] = 0; massU[i][j] = 0;}
        }
    massL[0][0] = massA[0][0];
    for (int i = 1; i < SLAR; ++i)
    {
        massU[0][i] = massA[0][i] / massL[0][0];
        massL[i][0] = massA[i][0];
    }
    for (int i = 1; i < SLAR; ++i)
    {
        massL[i][0] = massA[i][0];
        for (int k = 1; k < SLAR; ++k)
        {
            if (k <= i)
            {
                for (int m = 0; m < k; ++m)
                    massL[i][k] += (massL[i][m]*massU[m][k]);
                massL[i][k] = massA[i][k] - massL[i][k];
            }
            else
            {
                for (int m = 0; m < i; ++m)
                    massU[i][k] += (massL[i][m]*massU[m][k]);
            massU[i][k] = (massA[i][k] - massU[i][k]) / massL[i][i];
            }
        }
    }
    printf ("Press 1 if you want to see step output of LU-method\n");
    char ch = _getch();
    if (ch == '1') PrintLU (SLAR, massL, massU);
    for (int i = 0; i < SLAR; ++i) //визначаємо massY
    {
        double x = 0;
        for (int j = 0; j < i; ++j)
            x += massL[i][j]*massX[j];
        massX[i] = (massB[i] - x) / massL[i][i];
    }
    for (int i = SLAR-1; i >= 0; --i) //визначаємо корені рівняння
    {
        double x = 0;
        for (int j = SLAR-1; j > i; --j)
            x += massU[i][j]*massX[j];
        massX[i] = massX[i] - x;
    }
    PrintX(SLAR, massX);
    for (int i = 0; i < SLAR; ++i)
        free (massL[i]);
    free (massL);
    massL = NULL;
    for (int i = 0; i < SLAR; ++i)
        free (massU[i]);
    free (massU);
    massU = NULL;
    free (massX);
    massX = NULL;
}

void Jacobi(const int SLAR, double **massA, double *massB, double EPS)
{
    double *mass;
    CreateMass1(SLAR, &mass);
    double *massX;
    CreateMass1(SLAR, &massX);
    double *mass2;
    CreateMass1(SLAR, &mass2);
    double **mass1 = CreateMass2(SLAR, mass1);
    for (int i = 0; i < SLAR; ++i)
    {
        mass[i] = massB[i];
        mass2[i] = massB[i];
        for (int j = 0; j < SLAR; ++j)
        {
            mass1[i][j] = massA[i][j];
        }
    }
    bool check_on_eps; //змінна для перевірки точності
    printf ("Press 1 if you want to see step output of method of Jacobi\n\n");
    char ch = _getch();
    if (ch == '1')
    {
        for (int i = 0; i < SLAR; ++i) //виводимо систему на екран
        {
            for (int j = 0; j < SLAR; ++j)
                printf ("%+.2lfx[%i]", mass1[i][j], j+1);
        printf ("= %+.2lf\n", mass[i]);
        }
        putchar ('\n');
    }
    for (int i = 0; i < SLAR; ++i) //перетворюємо систему  до вигляду Х = ах+b
    {
        for (int j= 0; j < SLAR; ++j)
        {
            if (i == j) continue;
            mass1[i][j] = -mass1[i][j] / mass1[i][i];
        }
        mass[i] = mass2[i] / mass1[i][i];
        mass2[i] = mass2[i] / mass1[i][i];
        mass1[i][i] = 0;
    }
    if (ch == '1')
    {
        for (int i = 0; i < SLAR; ++i)
        {
            printf ("x[%i] = ", i+1);
            for (int j = 0; j < SLAR; ++j)
                if (j != i) printf ("%+.2lfx[%i]", mass1[i][j], j+1);
            printf ("%+.3lf\n", mass[i]);
        }
        putchar ('\n');
    }
    int k = 0; //змінна для кількості ітерацій
    do //метод Якобі
    {
        k++;
        for (int i = 0; i < SLAR; ++i)//дублюємо масиви, з яких mass лишається незмінним, а massX переприсваюємо
        {
            massX[i] = mass[i];
            mass[i] = mass2[i]; //для знаходження х треба додати вільний член mass2, тому одразу це врахуємо
        }
        for (int i = 0; i < SLAR; ++i) //знаходження нових значень для х
            for (int j = 0; j < SLAR; ++j)
                mass[i] += mass1[i][j]*massX[j]; //особливість методу Зейгеля від Якобі
        if (ch == '1') //виводимо проміжні значення на певній ітерації
        {
            printf ("Number of iteration for method of Jacobi: %d\n", k);
            for (int i = 0; i < SLAR; ++i)
                printf ("x[%d] = %+.10lf\n", i+1, mass[i]);
            putchar ('\n');
        }
        check_on_eps = true;
        for (int i = 0; i < SLAR; ++i) //перевіряємо досяжність заданої точності
            if (fabs(mass[i]-massX[i]) > EPS)
                check_on_eps = false;
    } while (check_on_eps == false);
    if (ch != '1')
    {
        putchar ('\n');
        for (int i = 0; i < SLAR; ++i)
            printf ("x[%d] = %+.10lf\n", i+1, mass[i]);
        putchar ('\n');
    }
    free (mass);
    mass = NULL;
    free (massX);
    massX = NULL;
    for (int i = 0; i < SLAR; ++i)
        free (mass1[i]);
    free (mass1);
    free (mass2);
}

void Zeidel(const int SLAR, double **massA, double *massB, double EPS)
{
    double *mass;
    CreateMass1(SLAR, &mass);
    double *massX;
    CreateMass1(SLAR, &massX);
    double *mass2;
    CreateMass1(SLAR, &mass2);
    double **mass1 = CreateMass2(SLAR, mass1);
    for (int i = 0; i < SLAR; ++i)
    {
        mass[i] = massB[i];
        mass2[i] = massB[i];
        for (int j = 0; j < SLAR; ++j)
        {
            mass1[i][j] = massA[i][j];
        }
    }
    bool check_on_eps; //змінна для перевірки точності
    printf ("Press 1 if you want to see step output of method of Zeidel\n\n");
    char ch = _getch();
    if (ch == '1')
    {
        for (int i = 0; i < SLAR; ++i) //виводимо систему на екран
        {
            for (int j = 0; j < SLAR; ++j)
                printf ("%+.2lfx[%i]", mass1[i][j], j+1);
        printf ("= %+.2lf\n", mass[i]);
        }
        putchar ('\n');
    }
    for (int i = 0; i < SLAR; ++i) //перетворюємо систему  до вигляду Х = ах+b
    {
        for (int j= 0; j < SLAR; ++j)
        {
            if (i == j) continue;
            mass1[i][j] = -mass1[i][j] / mass1[i][i];
        }
        mass[i] = mass2[i] / mass1[i][i];
        mass2[i] = mass2[i] / mass1[i][i];
        mass1[i][i] = 0;
    }
    if (ch == '1')
    {
        for (int i = 0; i < SLAR; ++i)
        {
            printf ("x[%i] = ", i+1);
            for (int j = 0; j < SLAR; ++j)
                if (j != i) printf ("%+.2lfx[%i]", mass1[i][j], j+1);
            printf ("%+.3lf\n", mass[i]);
        }
        putchar ('\n');
    }
    int k = 0; //змінна для кількості ітерацій
    do //метод Зейделя
    {
        k++;
        for (int i = 0; i < SLAR; ++i)//дублюємо масиви, з яких mass лишається незмінним, а massX переприсваюємо
        {
            massX[i] = mass[i];
            mass[i] = mass2[i]; //для знаходження х треба додати вільний член mass2, тому одразу це врахуємо
        }
        for (int i = 0; i < SLAR; ++i) //знаходження нових значень для х
            for (int j = 0; j < SLAR; ++j)
                if (j > i)  mass[i] += mass1[i][j]*massX[j]; //особливість методу Зейделя від Якобі
                else mass[i] += mass1[i][j]*mass[j];         //використання  (n+1)-го наближення
        if (ch == '1') //виводимо проміжні значення на певній ітерації
        {
            printf ("Number of iteration for method of Zeidel: %d\n", k);
            for (int i = 0; i < SLAR; ++i)
                printf ("x[%d] = %+.10lf\n", i+i, mass[i]);
            putchar ('\n');
        }
        check_on_eps = true;
        for (int i = 0; i < SLAR; ++i) //перевіряємо досяжність заданої точності
            if (fabs(mass[i]-massX[i]) > EPS)
                check_on_eps = false;
    } while (check_on_eps == false);
    if (ch != '1')
    {
        putchar ('\n');
        for (int i = 0; i < SLAR; ++i)
            printf ("x[%d] = %+.10lf\n", i+1, mass[i]);
        putchar ('\n');
    }
    free (mass);
    mass = NULL;
    free (massX);
    massX = NULL;
    for (int i = 0; i < SLAR; ++i)
        free (mass1[i]);
    free (mass1);
    free (mass2);
}

double InterpolateLagrangePolynomial (const double x, double* X, double* Y, int size)
{
    double y = 0; //значення, яке приймає функція у точці х
    for (int i = 0; i < size; i++)
    {
        double y1 = 1; //додаткова змінна для підрахунку окремих доданків
        for (int j = 0; j < size; j++)
        {
            if (j == i) continue;
            y1 *= (x - X[j])/(X[i] - X[j]);
        }
        y += y1*Y[i];
    }
    printf ("y = %lf\n", y);
    return y;
}

double InterpolateNewtonPolynomial(const double x, double* X, double* Y, const int size)
{
    double y = Y[0]; //задаємо шуканому значенню значення в першій точці
    double y1, y2;
    for (int i = 1; i < size; i++)
    {
        y1 = 0;
        for (int j = 0; j <= i; j++) //шукаємо f(x0,x1,...,xi)
        {
            y2 = 1;
            for (int k = 0; k <= i; k++)
            {
                if (k != j) y2 *= (X[j] - X[k]);
            }
            y1 += Y[j] / y2;
        }
        for (int k = 0; k < i; k++) y1 *= (x - X[k]);//шукаємо добутки (х-хі)
        y += y1;
    }
    printf ("y = %lf\n", y);
    return y;
}

void MinSquare(int iPoints, double* X, double* Y, int iPolinom)
{
    double **massA = new double*[iPolinom+1];
    double *massB = new double[iPolinom+1];
    if (!massA) NotEnoughMemory();
    if (!massB) NotEnoughMemory();
    for (int i = 0; i < iPolinom+1; ++i)
    {
        massA[i] = new double[iPolinom+1];
        if (!massA[i]) NotEnoughMemory();
    }
    massA[0][0] = iPoints;
    for (int i = 1; i < iPolinom+1; ++i) //заповнюємо верхню ліву половину матриці
    {
        double temp = 0; //змінна для підрахунку суми кожного елемента матриці
        for (int j = 0; j < iPoints; ++j)
        {
            temp += pow(X[j], (double) i);
        }
        int col = i, row = 0;
        do //заповнюємо всі діагональні елементи
        {
            massA[row][col] = temp;
            ++row; --col;
        } while (col != -1);
    }

    for (int i = 1; i < iPolinom+1; ++i) //заповнюємо нижню праву половину матриці
    {
        double temp = 0; //змінна для підрахунку суми кожного елемента матриці
        for (int j = 0; j < iPoints; ++j)
        {
            temp += pow(X[j], (double) i+iPolinom);
        }
        int col = i, row = iPolinom;
        do //заповнюємо всі діагональні елементи
        {
            massA[row][col] = temp;
            --row; ++col;
        } while (col != iPolinom+1);
    }

    for (int i = 0; i < iPolinom+1; ++i) //цикл для визначення значень вільних членів
    {
        massB[i] = 0;
        for (int j = 0; j < iPoints; ++j)
        {
            massB[i] += Y[j]*pow(X[j], i);
        }
    }
    Gauss(iPolinom+1, massA, massB);
    for (int i = 0; i < iPolinom+1; ++i)
    {
        delete (massA[i]);
    }
    delete (massA);
    massA = NULL;
    delete (massB);
    massB = NULL;
}

void Euler(const double x0, const double y0, const double x1, const double h)
{
    double x = x0, y = y0;
    cout << "   Euler's method\nx = " << x << "\t y = " << y << endl;
    while (x <= x1)
    {
        double f = function(x, y); //визначення приросту
        y += h*f;
        x += h;
        cout << "x = " << x << "\t y = " << y << endl;
    }
}

void Runge_Kutta(const double x0, const double y0, const double x1, const double h)
{
    double x = x0, y = y0;
    cout << "\nRunge-Kutta's method\nx = " << x << "\t y = " << y << endl;
    while (x <= x1)
    {
        double k1 = h* function(x,y); //визначення коефіцієнтів
        double k2 = h * function(x+h/2, y+k1/2);
        double k3 = h * function(x+h/2, y+k2/2);
        double k4 = h * function(x+h, y+k3);
        double f = (k1+2*k2+2*k3+k4)/6; //визначення приросту
        y += f;
        x += h;
        cout << "x = " << x << "\t y = " << y << endl;
    }
}

void Adams(const double x0, const double y0, const double x, const double h)
{
    double x4 = x0, y = y0, y4 = y0, x1 = x4+3*h, x2 = x4 + 2 * h, x3 = x4 + h;
    double y3 = y4 + h * function(x4, y4);
    double y2 = y3 + h*(3*function(x3, y3) - function(x4, y4))/2;
    double y1 = y2 + h*(23*function(x2,y2)-16*function(x3, y3) + 5*(x4, y4))/12;
    cout << "\n   Adams method\nx = " << x4 << "\t y = " << y4 << endl;
    cout << "x = " << x3 << "\t y = " << y3 << endl;
    cout << "x = " << x2 << "\t y = " << y2 << endl;
    cout << "x = " << x1 << "\t y = " << y1 << endl;
    while (x1 <= x)
    {
        y = y1 + h*(55*function(x1,y1)-59*function(x2,y2)+37*function(x3,y3)-9*function(x4,y4))/24;
        y4 = y3;
        y3 = y2;
        y2 = y1;
        y1 = y;
        x4 = x3;
        x3 = x2;
        x2 = x1;
        x1 += h;
        cout << "x = " << x1 << "\t y = " << y << endl;
    }
}