/*Для тестирования C#-программы принят пример 2.3, приведенный в книге Шуп. Т. 
Решение инженерных задач на ЭВМ: Практическое ру-ководство: 
Пер. с англ. – М.: Мир, 1982. – С.41-43. */
using System; //C#-программа
using System.IO;
/* Метод Ньютона решения системы НАУ - 
Method Newton decisions of the SNAE (four equations).
Decision of the system of the nonlinear algebraic equations (NAE)*/

namespace ConsAppl_SNAE_4
{
    class DecisionSNAE
    {
        public DecisionSNAE() { }

        // Реализуется метод Гаусса
        private double[] GAUS(double[,] A, double[] B, double[] X, int n)
        {
            double factor, temp, adding;
            int kp1, ip1, nm1, l, i;

            nm1 = n - 2;
            for (int k = 0; k <= nm1 + 1; k++)
            {
                kp1 = k + 1;
                //перестановка строк
                l = k;

                for (i = kp1; i < n; i++)
                {
                    if ((Math.Abs(A[i, k]) - Math.Abs(A[l, k])) > 0) l = i;
                }

                if ((l - k) > 0)
                {
                    for (int j = k; j < n; j++)
                    {
                        temp = A[k, j]; A[k, j] = A[l, j]; A[l, j] = temp;
                    }
                    temp = B[k]; B[k] = B[l]; B[l] = temp;
                }

                {
                    // исключение неизвестных - приведение системы к треугольному виду; прямой ход

                    for (i = kp1; i < n; i++)
                    {
                        factor = A[i, k] / A[k, k]; A[i, k] = 0;
                        for (int j = kp1; j < n; j++)
                        {
                            A[i, j] = A[i, j] - factor * A[k, j];
                        }
                        B[i] = B[i] - factor * B[k];
                    }
                }
            }

            // определение неизвестных; обратный ход
            X[n - 1] = B[n - 1] / A[n - 1, n - 1]; i = nm1;
            while (i >= 0)
            {
                ip1 = i + 1; adding = 0;
                for (int j = ip1; j < n; j++)
                {
                    adding = adding + A[i, j] * X[j];
                }
                X[i] = (B[i] - adding) / A[i, i]; i = i - 1;
            }
            return X;
        }

        private static void Main()
        {
            int n = 4; int itmax = 50; double eps = 0.0001;

            double[] X = new double[4] { 1.0, 1.0, 1.0, 1.0 };
            double[] fv = new double[4] { 0.0, 0.0, 0.0, 0.0 };
            double[] dx = new double[4] { 0.0, 0.0, 0.0, 0.0 };
            double[,] fd = new double[4, 4] { { 0.0, 0.0, 0.0, 0.0 },
                                                              { 0.0, 0.0, 0.0, 0.0 },
                                                              { 0.0, 0.0, 0.0, 0.0 },
                                                              { 0.0, 0.0, 0.0, 0.0 } };

            DecisionSNAE metodN = new DecisionSNAE();

            FileStream file = new FileStream("result.txt", FileMode.Append);
            StreamWriter write = new StreamWriter(file);

            write.WriteLine();
            write.WriteLine("Решение системы {0}-x НАУ методом Ньютона.", n);
            write.WriteLine();
            write.WriteLine("Заданная точность решения            eps = {0}", eps);
            write.WriteLine("Максимальное число итераций       itmax = {0}", itmax);
            write.WriteLine("Начальные значения корней: ");
            for (int j = 0; j < n; j++) write.Write("X({0})={1:#.##}\t", j + 1, X[j]); write.WriteLine();
            write.WriteLine();
            write.WriteLine("Результаты решения системы НАУ: ");

            for (int i = 0; i < itmax; i++)
            {// метод Ньютона
                metodN.Ffv(X, fv);
                metodN.Ffd(X, fd);
                for (int j = 0; j < n; j++) if (Math.Abs(fv[j]) > eps) goto loop2;
                goto loop1;
                loop2: metodN.GAUS(fd, fv, dx, n);
                for (int k = 0; k < n; k++) X[k] = X[k] + dx[k]; write.Write("i ={0}  ", i + 1);
                for (int j = 0; j < n; j++) write.Write("\tX({0})={1:#.####}", j + 1, X[j]); write.WriteLine();
                loop1:;
            }
            write.Flush();
            write.Close();
            file.Close();
        }
        private double[,] Ffd(double[] X, double[,] fd)
        {
            // элементы якобиана 
            fd[0, 0] = 1; fd[0, 1] = 2;
            fd[0, 2] = 1; fd[0, 3] = 4;
            fd[1, 0] = 2 * (X[0] + X[1]); fd[1, 1] = 2 * X[0];
            fd[1, 2] = 0; fd[1, 3] = 3 * X[3] * X[3];
            fd[2, 0] = 3 * X[0] * X[0]; fd[2, 1] = 0;
            fd[2, 2] = 2 * X[2]; fd[2, 3] = 1;
            fd[3, 0] = 0; fd[3, 1] = 3;
            fd[3, 2] = X[3]; fd[3, 3] = X[2];
            return fd;
        }
        private double[] Ffv(double[] X, double[] fv)
        {  //элементы вектора невязки
            fv[0] = -X[0] - 2 * X[1] - X[2] - 4 * X[3] + 20.7;
            fv[1] = -X[0] * X[0] - 2 * X[0] * X[1] - X[3] * X[3] * X[3] + 15.88;
            fv[2] = -X[0] * X[0] * X[0] - X[2] * X[2] - X[3] + 21.218;
            fv[3] = -3 * X[1] - X[2] * X[3] + 21.1;
            return fv;
        }
    }
}

