0

Всем добрый день! Стоит задача решить слау, при этом число строк и столбцов не обязательно одинаковое (решаю задачу с тимуса 1580: http://acm.timus.ru/problem.aspx?space=1&num=1580). Матрица коэффициентов имеет значения только 1 или 0, а правые части могут быть любыми. Это завело меня в тупик связанный с тем, что на некоторые входные данные алгоритм дает NaN. Я разобрался что это из-за того что при приведении к треугольному виду используется деление на 0. Подскажите пожалуйста, как обойти это? Я нашел что есть т.н. метод Барейса, но литературы про него крайне мало. Помогите исправить код, буду очень благодарен:

using System;
class Program
{
    static void Main(string[] args)
    {
        //INPUT
        string[] inLine = Console.ReadLine().Split();
        int N = int.Parse(inLine[0]);   //N - число неизвестных (число студентов)
        int M = int.Parse(inLine[1]);   //M - число уравнений (число пар студентов пришедших к декану)
        double[,] A = new double[M, N+1];   //A - матрица коэффициентов
        for (int i = 0; i < M; i++)
        {
            inLine = Console.ReadLine().Split();
            A[i, int.Parse(inLine[0]) - 1] = 1;
            A[i, int.Parse(inLine[1]) - 1] = 1;
            A[i, N] = int.Parse(inLine[2]);
        }
        if (M < N)
            Console.WriteLine("IMPOSSIBLE");
        else
        {
            double[] result = GaussMethod(A);
            for (int i = 0; i < M; i++)
                if (double.IsNaN(result[i]))
                {
                    Console.WriteLine("IMPOSSIBLE");
                    return;
                }
            for (int j = 0; j < result.Length; j++) 
                Console.WriteLine("{0:0.00}", result[j]);
        }
    }

    static double[] GaussMethod(double[,] a)
    {
        //найти размер матрицы
        int width = a.GetLength(1);
        int height = a.GetLength(0);

        double[] x = new double[height]; //массив решений уравнения
        //инициализация массива корней уравнения
        for (int i = 0; i < x.Length; i++) 
        {
            x[i] = a[i, width - 1];
        }

        //приведение матрицы к треугольному виду
        double m;
        for (int i = 0; i < width - 1; i++) 
            for (int j = i + 1; j < height; j++) 
            {
                m = a[j, i] / a[i, i];
                for (int k = 0; k < height; k++) 
                    a[j, k] -= m * a[i, k]; 
                x[j] -= m * x[i];
            }

        for (int i = 0; i < height; i++) 
            a[i, width - 1] = x[i];


        for (int i = height - 1; i >= 0; i--) 
        {
            for (int j = i + 1; j < width - 1; j++) 
            {
                x[i] -= i != j ? a[i, j] * x[j] : 0;
            }
            x[i] = x[i] / a[i, i];
        }           
        return x;
    }
}
3
  • а чем вам метод Гаусса не угодил? Реализуется легко, только надо не забывать переставлять уравнения в матрице. Тогда и деления на ноль не будет.
    – rdorn
    18 апр 2016 в 19:26
  • а как именно их переставлять? если видим что грядет ноль в знаменателе, то ищем ту где нуля в этом столбце нет и меняем? а что если такой строки нет? тогда не совместна система?
    – I.R.
    18 апр 2016 в 20:54
  • ну вот, теорию вроде знаете, осталось сделать =) хотя, по поводу не совместности системы, не уверен, уточните.
    – rdorn
    18 апр 2016 в 21:00

2 ответа 2

2

Внимательнее посмотрите на вот этот участок кода:

    for (int i = 0; i < width - 1; i++) 
        for (int j = i + 1; j < height; j++) 
        {
            m = a[j, i] / a[i, i];

Вы не проверяете знаменатель на 0 - отсюда и NaNы. Вам нужно среди всех строк в качестве опорной выбрать ту, которая не имеет ноля в i-м столбце, после чего поменять ее с i-й строкой.

1

Очевидно, почему Ваш алгоритм даёт на некоторые данные значения NAN. В теории численных решений СЛАУ есть такое понятие, как устойчивость системы. Вот именно здесь возникает Ваша проблема судя по всему. Для начала, я расскажу это на примере LU-разложения и LUP-разложения. Я полагаю, что Ваш метод Барейса здесь ни при чём! Как я вижу, в вашем решении реализуется именно LU-разложение. А нужно Вам LUP-разложение. Попробуйте реализовать его. Если не зайдёт, пишите.

Перед тем как начать, скажем, что существуют такие системы, в которых значения в матрице сильно различаются. Например:

1    0    1
1000 1    1000
1    1000 500

Разумеется, с аналитической точки зрения всё у нас прекрасно. Но с точки зрения вычислений -- всё плохо. Если рассматривать учитывать тот факт, что накапливается погрешность.

Понятно, что от перемены мест уравнений ничего не изменится с аналитической точки зрения. Но от порядка операций зависит очень много, коль скоро мы говорим о нарастании погрешности. Поэтому при численном решении СЛАУ на машине порядок следования уравнений очень важен!

Теперь, для того, чтобы продемонстрировать, что такое устойчивость, рассмотрим 2 уравнения (здесь Вас должен постигнуть катарсис):

Неустойчивая СЛАУ Фиг. 1

Тут же рассмотрим вторую СЛАУ:

Устойчивая СЛАУ Фиг. 2

Видно, что чем ближе к параллельности прямые, тем меньше будут отличаться коэффициенты. То же самое будет и в пространствах большей размерности. Другое дело, что это представить нельзя для n > 3.

Задачу об определении "качества" матрицы также формулируют, как вопрос о плохой обусловленности матрицы. Это понятие тесно связано с понятием нормы матрицы в некотором нормированном пространстве. Норма -- это обобщение понятия длины вектора для абстрактных пространств (есть также, схожее определение, называемое метрикой -- длина). Вообще говоря, в многомерных пространствах понятие вектора и точки смазываются. Это, вероятно, происходит потому что точка -- это неудобный объект. Его сложно обобщить. Вектор же обобщить проще. Поэтому часто говорят о точке, как о радиус-векторе. В таком случае, можно искать и расстояние между векторами при помощи метрики. Таким образом, метрика -- расстояние. Норма -- длина.

Разные нормы используются для разных вещей. Например, max-норма (inf-норма) используется для анализа выбросов в решении. Евклидова норма хорошо подходит для анализа решения в целом (так сказать, среднее по больнице). Она плохо реагирует на выбросы.

Для того, чтобы понять, насколько матрица хорошо или плоха, вводят так называемые числа обусловленности, которые можно точно посчитать аналитически. Число обусловленности. Но в силу того, что их сложно считать (вычисление обратной матрицы очень трудоёмкая операция), часто прибегают к вычислению их оценок.

Когда я разбирался с этой темой, у меня возникал вопрос: "Почему обусловленность матрицы считается именно так?" Первый ответ -- оно хорошо характеризует величину "хорошести" матрицы. Для второго ответа нужно сказать, что часто на практике бывает так, что у одной и то же системы бывают разные правые части, т.е.

Ax = b

и

Ax* = c

В таком случае, возникает вопрос, насколько сильно будет отличаться СЛАУ, если

c = b + epsilon,

где epsilon -- бесконечно малая величина. Вот здесь и выплывает второй ответ. Берут отношение (здесь тоже Вас тоже должен постигнуть катарсис)

x*/x,

которое и характеризует "качество" матрицы. Если Вы хотите подробнее разобраться с этой темой, можете посмотреть в сторону теории регуляризации и теории устойчивости по Ляпунову.

Вернёмся к нашим баранам. Основным косяком, который есть в LU -- это деление очень маленького числа на очень большое. Ведь именно в результате такой операции возникают ещё меньшие числа, чем были в начале, которые уже могут чпокнуться машиной.

Чтобы этого не происходило, нужно сделать очевидную вещь -- поставить самые большие числа в числитель. Таким образом, мы получим куда лучшее решение, чем было до этого. Не так ли!? Но будьте внимательны! Вам нужно запоминать, какие вектора вы переставляли, так как при перестановке уравнений, меняется и полученный в ответе вектор! Поэтому нужно будет переставить всё в обратном порядке.

3
  • У автора нет проблем с "почти параллельными" прямыми - у него же только 1 и 0 в матрице... У автора проблемы с делением на ноль же! 1 июн 2016 в 5:40
  • Деление на ноль (почти на ноль) -- всё та же проблема обусловленности, не? 1 июн 2016 в 6:14
  • Это частный случай, который можно объяснить проще 1 июн 2016 в 7:05

Ваш ответ

By clicking “Отправить ответ”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Всё ещё ищете ответ? Посмотрите другие вопросы с метками или задайте свой вопрос.