using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;

/* ЗАДАНИЕ
F(t,X) = { X1 + Exp(t), X1*X1 + 2*X2, Sqrt(X3), X4, X4*X5 }
X0 = { 1, 1, 1, 1, 1 }
t0 = 0;
tk = 1;
*/

namespace ConsoleApplication22
{
    class Program
    {
        const int n = 5;

        static double[,] U = new double[n, n];
        static double[,] L = new double[n, n];
        static double[,] A = new double[n, n];
        static double[] B = new double[n];
        static double[] Z = new double[n];

        static double[] X = new double[n];
        static double[,] M = new double[n, n];
        static double[] dX = new double[n];
        static double[] Xo = { 1, 1, 1, 1, 1 };
        static double[] f = new double[n];
        static double[] F = new double[n];
        static int i, j, k, m;
        static double s, to = 0, tk = 1, dt = 0.1, t;
        static void Main(string[] args)
        {
            while (t < tk)
            {
                t = to + dt;
                fun1();
                for (i = 0; i < n; ++i)
                {
                    X[i] = Xo[i] + dt * f[i];
                }
                function_lab2();
                for (i = 0; i < n; ++i)
                {
                    Xo[i] = X[i];
                }
                Console.WriteLine("X1 = " + X[0] + "; " + "X2 = " + X[1] + "; " + "X3 = " + X[2] + "; " + "X4 = " + X[3] + ";"); //X[1], X[2], X[3], X[4]);
                to = t;
            }
            Console.Read();
        }

        static void fun1()
        {
            f[0] = Xo[0] + Math.Exp(t);
            f[1] = Xo[0] * Xo[0] + 2 - Xo[1];
            f[2] = Math.Sqrt(Xo[2]);
            f[3] = Xo[3];
            f[4] = Xo[3] * Xo[4];
        }

        static void function()
        {
            F[0] = X[0] - Xo[0] - dt * f[0];
            F[1] = X[1] - Xo[1] - dt * f[1];
            F[2] = X[2] - Xo[2] - dt * f[2];
            F[3] = X[3] - Xo[3] - dt * f[3];
            F[4] = X[4] - Xo[4] - dt * f[4];
        }

        static void FindLU()
        {
            for (i = 0; i < n; i++)
            {
                for (j = 0; j < n; j++)
                {
                    U[i,j] = 0;
                    L[i,j] = 0;
                }
            }
            for (i = 0; i < n; i++)
            {
                U[i,i] = 1;
            }
            k = 0;

            for (i = 0; i < n; i++)
            {
                L[i,0] = A[i,0];
            }

            for (j = 0; j < n; j++)
            {
                U[0,j] = A[0,j] / L[0,0];
            }
            while (k < n - 2)
            {
                k++;
                for (i = k; i < n; i++)
                {
                    s = 0;
                    for (m = 0; m < k; m++)
                    {
                        s = s + L[i,m] * U[m,k];
                    }
                    L[i,k] = A[i,k] - s;
                }
                for (j = (k + 1); j < n; j++)
                {
                    s = 0;
                    for (m = 0; m < k; m++)
                    {
                        s = s + L[k,m] * U[m,j];
                    }
                    U[k,j] = (A[k,j] - s) / L[k,k];
                }
            }
            s = 0;
            for (m = 0; m < (n - 1); m++)
            {
                s = s + L[n - 1, m] * U[m, n - 1];
            }
            L[n - 1, n - 1] = A[n - 1, n - 1] - s;

        }

        static void FindZX()
        {
            Z[0] = B[0] / L[0,0];

            for (k = 1; k < n; k++)
            {
                s = 0;
                for (j = 0; j < k; j++)
                {
                    s = s + L[k,j] * Z[j];
                }
                Z[k] = (B[k] - s) / L[k,k];
            }
            dX[n - 1] = Z[n - 1];
            for (k = (n - 2); k > -1; k--)
            {
                s = 0;
                for (j = (k + 1); j < n; j++)
                {
                    s = s + U[k,j] * dX[j];
                }
                dX[k] = Z[k] - s;
            }
        }

        static void matrixJakobi()
        {
            M[0,0] = 1 - dt;
            M[0,1] = 0;
            M[0,2] = 0;
            M[0,3] = 0;
            M[0,4] = 0;

            M[1,0] = 2 * X[0] * dt;
            M[1,1] = 1 + dt;
            M[1,2] = 0;
            M[1,3] = 0;
            M[1,4] = 0;

            M[2,0] = 0;
            M[2,1] = 0;
            M[2,2] = 1 - (1 * dt) / (2 * Math.Sqrt(X[2]));
            M[2,3] = 0;
            M[2,4] = 0;

            M[3,0] = 0;
            M[3,1] = 0;
            M[3,2] = 0;
            M[3,3] = 1 - dt;
            M[3,4] = 0;

            M[4,0] = 0;
            M[4,1] = 0;
            M[4,2] = 0;
            M[4,3] = -dt;
            M[4,4] = 1 + dt;
        }

        static void function_lab2()
        {
            int stop = 0;
            while (true)
            {
                function();
                stop = 1;

                for (i = 0; i < n; i++)
                {
                    if (Math.Abs(F[i]) > 0.0001)
                    {
                        stop = 0;
                    }
                }

                if (stop != 0) break;
                matrixJakobi();
                for (i = 0; i < n; i++)
                {
                    B[i] = -F[i];
                }

                for (i = 0; i < n; i++)
                {
                    for (j = 0; j < n; j++)
                    {
                        A[i,j] = M[i,j];
                    }
                }
                FindLU();
                FindZX();

                for (i = 0; i < n; i++)
                {
                    X[i] = X[i] + dX[i];
                }
            }

        }
    }
}