1

Тема: Степеневий алгоритм.РМ-алгоритм

Добрий вечір. Реалізовую РМ-алгоритм, і роблячи перевірку результат не сходиться , не можу зрозуміти в чому проблема
https://imgur.com/7Ychles
https://imgur.com/OLq1n2V
Перевірка:Ax=hx
Початкові дані:
Матриця:2 1 0|1 3 1|0 1 2
h0=y0=1 1 1
eps=b=0.001
Перевірка може відбуватися з невеликою похибкою але в моєму результаті вона надто велика
Прошу допомогти, можливо неправильно рахую, або десь є якась помилка, буду дуже вдячний за допомогу

Прихований текст
class Program
    {
        public static double[] mat_vect(double[,] mat, double[] vect)
        {
            double[] check = new double[vect.Length];
            for (int h = 0; h < vect.Length; h++)
            {

                check[h] = 0;
                for (int k = 0; k < vect.Length; k++)
                    check[h] += (mat[h, k] * vect[k]);
            }
            return check;

        }
        public static double norm_y(double[] y)
        {
            double sum = 0;
            double norm = 0;
            for(int i = 0; i < y.Length; i++)
            {
                sum += Math.Pow(y[i], 2);
            }
            norm = Math.Sqrt(sum);
            return norm;
        }
        public static double norm(double[] x0, double[] x1)
        {
            double max = 0;
            List<double> m = new List<double>();
            for (int i = 0; i < x0.Length; i++)
            {
                m.Add(Math.Abs(x1[i] - x0[i]));
            }
            max = m.Max();
            return max;
        }
        public static void pm(double[,] mat, double[] h0, double[] y0, double eps, double b, int m)
        {
            int k = 0;
            double norma = norm_y(y0);//норма у0
            double[] x0 = new double[m];
            double[] h = new double[m];
            double[] x = new double[m];
            double[] y = new double[m];
            List<int> s0 = new List<int>();//множина індексів
            for (int i = 0; i < m; i++)
            {
                x[i] = (y0[i] / norma);//пункт 2 алгоритму, вектор х0
            }
            List<int> s = new List<int>();
            for (int i = 0; i < m; i++)
            {
                if (Math.Abs(x[i]) > b)
                {
                    s.Add(i);
                }
            }
            do
            {
                s0 = s;
                for(int i = 0; i < m; i++)
                {
                    x0[i] = x[i];
                }
                k++;
                y = mat_vect(mat, x0);//пункт 3, обчислення вектору у
                for (int i = 0; i < m; i++)
                {
                    h[s0[i]] = (y[s0[i]] / x0[s0[i]]);//пункт 4, вектор h
                }
                for (int i = 0; i < m; i++)
                {
                    x[i] = (y[i] / norm_y(y));//пункт 5, новий вектор х
                }
                s.Clear();
                for (int i = 0; i < m; i++)
                {
                    if (Math.Abs(x[i]) > b)
                    {
                        s.Add(i);
                    }
                }
            }
            while ((norm(h, h0) > eps) && (norm(x, x0) > eps));//пункт 6, перевірка критерію зупинки
            double sum = 0;
            for (int i = 0; i < h.Length; i++)
            {
                sum += h[i];
            }
            double sk = 0;
            for (int i = 0; i < s.Count(); i++)
            {
                sk += Math.Pow(s[i], 2);
            }
            sk = Math.Sqrt(sk);
            double h1 = (sum / sk);
            Console.WriteLine("h="+h1);//виведення результатів
            Console.WriteLine("x:");
            for (int i = 0; i < m; i++)
            {
                Console.WriteLine(x[i]);
            }
            Console.WriteLine("k="+k);

            Console.WriteLine();
            Console.WriteLine();
            Console.WriteLine();
            Console.WriteLine("Ax:");//виведення перевірки для візуального порівняння
            double[] Ax = mat_vect(mat, x);
            for(int i = 0; i < m; i++)
            {
                Console.WriteLine(Ax[i]);
            }
            Console.WriteLine("hx:");
            for(int i = 0; i < m; i++)
            {
                double tmp = h1 * x[i];
                Console.WriteLine(tmp);
            }
        }
        static void Main(string[] args)
        {
            int n, m;
            double eps,b;
            Console.WriteLine("Type m:");
            m = Convert.ToInt32(Console.ReadLine());
            n = m;
            double[] h0 = new double[m];
            double[] y0 = new double[m];
            double[,] mat = new double[m, n];
            for (int i = 0; i < m; i++)
            {
                for (int j = 0; j < n; j++)
                {
                    Console.WriteLine($"M[{i + 1}][{j + 1}]=");
                    mat[i, j] = Convert.ToDouble(Console.ReadLine());
                }
            }

            for (int i = 0; i < m; i++)
            {
                Console.WriteLine($"h[{i + 1}]:");
                h0[i] = Convert.ToDouble(Console.ReadLine());//зчитування h0
            }
            for (int i = 0; i < m; i++)
            {
                Console.WriteLine($"y[{i + 1}]:");
                y0[i] = Convert.ToDouble(Console.ReadLine());//зчитування у0
            }
            Console.WriteLine("Type eps:");
            eps = Convert.ToDouble(Console.ReadLine());
            Console.WriteLine("Type b:");
            b = Convert.ToDouble(Console.ReadLine());
            pm(mat, h0, y0, eps, b, m);
            Console.ReadLine();
        }
    }