Точность вычислений. Метод Гаусса

[H g]

Есть у меня система линейных алгебраических уравнений

a11x1 + a12x2 + a13x3 + ... a1NxN = b1
a21x1 + a22x2 + a23x3 + ... a2NxN = b2
a31x1 + a32x2 + a33x3 + ... a3NxN = b3
...
aN1x1 + aN2x2 + aN3x3 + ... aNNxN = bN

Получаю марицы a[][], b[], x[]. Решаю с помощью метода Гаусса.
Все работает, но есть одна проблема: накапливается погрешность и в итоге вместо точного ответа в 10000, имею 9999,

Код привожу, но он возможно требует пояснений, так как вырван из контекста

Скрытый текст

void CalculateMatrixSolve()
{
    int i, j, k, n, p, r;
    float a[100][100];
    float b[100];
    float x[100];
    float temp, eps;
    n = e::calculatematrixubound; // число уравнений
    eps = 0.0000001f;
    for (i=0; i<e::calculatematrixshift; i++) // это по сути обертка - внешние исходные данные, типа расчет выполнить i раз с разными входными данными.
                                              // переменная i в самом алгоритме не используется
    {
        for (j=0; j<n; j++)
        {
            for (k=0; k<n; k++)
            {
                a[j][k] = e::calculatematrixa[i][j][k]; // входная матрица A[j][k]
            }
            b[j] = e::calculatematrixb[i][j]; // входная матрица B[j]
            x[j] = 0.0f;
        }
        for (j=0; j<n; j++)
        {
            if (a[j][j] == 1)
            {
                for (r=j+1; r<n; r++)
                {
                    if (fabs(a[r][j]) >= eps)
                    {
                        temp = a[r][j];
                        for (p=j; p<n; p++)
                        {
                            a[r][p] -= (a[j][p] * temp);
                        }
                        b[r] -= (b[j] * temp);
                    }
                }
            }
            else
            {
                if (fabs(a[j][j]) < eps)
                {
                    for (r=j+1; r<n; r++)
                    {
                        if (fabs(a[r][j]) >= eps)
                        {
                            for (p=j; p<n; p++)
                            {
                                temp = a[j][p];
                                a[j][p] = a[r][p];
                                a[r][p] = temp;
                            }
                            temp = b[j];
                            b[j] = b[r];
                            b[r] = temp;
                            break;
                        }
                    }
                    if (fabs(a[j][j]) >= eps)
                    {
                        temp = a[j][j];
                        for (p=j; p<n; p++)
                        {
                            a[j][p] /= temp;
                        }
                        b[j] /= temp;
                        for (r=j+1; r<n; r++)
                        {
                            if (a[r][j] != 0)
                            {
                                temp = a[r][j];
                                for (p=j; p<n; p++)
                                {
                                    a[r][p] -= (a[j][p] * temp);
                                }
                                b[r] -= (b[j] * temp);
                            }
                        }
                    }
                    else
                    {
                        return;
                    }
                }
                else
                {
                    temp = a[j][j];
                    for (p=j; p<n; p++)
                    {
                        a[j][p] /= temp;
                    }
                    b[j] /= temp;
                    for (r=j+1; r<n; r++)
                    {
                        if (a[r][j] != 0)
                        {
                            temp = a[r][j];
                            for (p=j; p<n; p++)
                            {
                                a[r][p] -= (a[j][p] * temp);
                            }
                            b[r] -= (b[j] * temp);
                        }
                    }
                }
            }
        }
        for (j=(n-1); j>=0; j--)
        {
            x[j] = b[j];
            for (p=0; p<j; p++)
            {
                b[p] -= (a[p][j] * x[j]);
            }
        }
        for (j=0; j<n; j++)
        {
            e::calculatematrixx[i][j] = x[j]; // выходная матрица X[j] результат
        }
    }
}

Здесь есть некоторая погрешность (eps), потому что сам метод предполагает деление, чтобы получилась "треугольная" матрица. Так вот, если делитель будет например 3, то 1.0/3.0 = 0.3333... значит деление не может быть выполнено точно и дальше накапливается погрешность. Попробовал поменять (eps) но было еще хуже, при чем при увеличении числа знаков после запятой функция решает неправильно. Путем практики подобрал (eps = 0.0000001f;)
// Если 1.0/3.0 полчим 0.333333, а если потом 0.333333 умножить на 3.0 то получим 0.999999 а не 1.0
// В этом мне кажется загвоздка

Вопрос как увеличить точность решения ?
Предполагаемые варианты:
1) заменить float на double, поможет ли это
2) использовать другой метод не Гаусса, тогда какой..

Может что-нибудь посоветуете ?

Сообщение отредактировано: H g - 07.06.24, 22:01

[Majestio]

Цитата H g @ 07.06.24, 21:35

Вопрос как увеличить точность решения ?

Интересный вопрос. Я когда-то им задавался, но не было прикладной задачи - и так руки и не дошли.

Тем не менее, поделюсь своими соображениями ...

Сперва немного рассуждений. Классический пример, приведенный выше (деление, а потом умножение на три), показывает, прямо тыкает носом в то, что использование float и double - тупиковый метод. Возникает законный вопрос, ачо делать-та? И тут возникает ответ: не работайте с десятичными дробями, работайте со смешанными дробями! Естественно, это нужно реализовывать самому (ну или поискать в сети, может уже кто удосужился).

Т.е. любое число у нас представляется примерно следующим образом:

class BigDecimal {
  string WholePart;   // целая часть
  string Numerator;   // числитель
  string Denominator; // знаменатель
}

Ну и далее реализация кучи публичных методов и перегрузок операций. Сразу возникает законный вопрос "WTF strings?". Рассказываю, может непоследовательно, но тезисно и с душой

1. В своей реализации вообще-то BigDecimal мы должны задействовать какую-либо реализацию BigInt. Я встречал простейшие реализации бигинтов на строках. Возможно эти части и следовало бы объявлять не строками, а сразу большими целыми, но не суть...
2. В процессе арифметических операций у нас вообще ничего не теряется ^*, Но за это нужно платить быстро растущими размерами всех полей BigDecimal, равно как и значительным снижением скорости вычислений
3. В результате нас практически перестает беспокоить потеря точности, но начинает серьёзно напрягать проблема потери производительности и потребления ресурсов
4. Поэтому, время от времени, по реализованной нами, вернее вами, стратегии, BigDecimal должен производить нормализацию своих данных путем поиска НОД числителя и знаменателя, и сокращать их на этот НОД
5. Ну и второй вариант - если нормализации явно недостаточно и данные все равно слишком большие, должна быть разработана стратегия управления уменьшением точности - сокращением числителя и знаменателя (одного из) делением на какое-то серьёзное целое нацело, а второго из ... с потерями. Но когда это касается, к примеру 30 000-й цифры после запятой, потери точности буду едва заметны. В любом случае для решения мелких и средних задач по вычислениям хватит и решения из п.4.
6. И конечно последний шаг - уже после всех вычислений мы можем спокойно вызвать метод, который из BigDecimal нам вернет float или double

^* с учетом п.5

Как резюме

В такой реализации мы сохраняли все что можно в процессе вычислений, платили мучительными секундами ожидания и повышения температуры центрального процессора, но получили на выходе максимально возможный по точности результат во float или double.

Перед реализацией своего BigDecimal нужно очень внимательно выбирать либу BigInt, т.к. реализация на строках последней - не лучший вариант. Чуйка подсказывает, что на битовых массивах будет получше.

Ну и последнее. Ложечка дегтя в вашу баночку с вашим мёдом. Ясный перец, что речь идет только об элементарных арифметических операциях. В BigDecimal мы не сможем хранить иррациональные числа, не сможем обеспечивать 100% точность при извлечении корней, вычислении логарифмов & etc. Жизнь - боль!

[H g]

Majestio Спасибо за подробный ответ.

Цитата

не работайте с десятичными дробями, работайте со смешанными дробями

В этом что-то есть, хорошая мысль. В целом если умножение/деление будет выполнено на последнем шаге, то все типа /3 и x3 в конце сократятся. Нужно подумать как реализовать.

[amk]

А я бы посоветовал для начала не сравнивать диагональный элемент с каким-то эпсилон, а всегда искать главный элемент. Достаточно искать по столбцу. Главный - это максимальный по величине (без учёта знака). Своим сравнением с эпсилон ты допускаешь, что коэффициент на который ты умножаешь ведущую строку при исключении переменной из остальных строк может оказываться много больше 1. А на этот коэффициент умножаются не только элементы матрицы, но и их погрешности.
Собственно алгоритмы без выбора главного элемента в программах давно уже не применяются. Они остались только на уроках математики при ручном решении СЛАУ, когда ученик глазами видит, что унего происходит.

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

Majestio, ты бы лучше посоветовал ему не изобретать велосипед, а прикрутить к программе какую-нибудь хорошо отлаженную библиотеку вычислений c повышенной точностью, gmp например.

[Majestio]

Цитата amk @ 08.06.24, 20:12

Majestio, ты бы лучше посоветовал ему не изобретать велосипед, а прикрутить к программе какую-нибудь хорошо отлаженную библиотеку вычислений c повышенной точностью, gmp например.

amk, ты бы лучше не мне советовал как советовать, а посоветовал бы свой совет ТС'у. А давай ты посоветуешь предложенную тобою либу? Заодно и покажешь ее качество и "невелосипедность" на примерах? А мы проверим.

Добавлено 09.06.24, 03:33
Хотя ... чего катать вату. Возьмем и накидаем пример деления и последующего умножения на "три", с использованием библиотеки GMP, очень отлаженной и очень надёжной. Ну и посмотрим качество совета.

#include <stdio.h>
#include <gmp.h>
 
int main(void) {
    mpz_t x, result;
    mpz_init_set_str(x, "10", 10);
    mpz_init(result);
 
    // Деление на 3
    mpz_tdiv_q_ui(x, x, 3);
 
    // Умножение на 3
    mpz_mul_ui(result, x, 3);
 
    gmp_printf("Результат: %Zd\n", result);
 
    mpz_clear(x);
    mpz_clear(result);
    return 0;
}

И так. Число 10 делим на 3, делим, делим, делим, а потом умножаем на то же 3. Барабаны, фейерверки ... и получаем число 9. Аплодисменты и занавес!

amk, а вывод простой - советовать другим как советовать - это явно не твое!

[amk]

Цитата H g @ 08.06.24, 06:27

В целом если умножение/деление будет выполнено на последнем шаге, то все типа /3 и x3 в конце сократятся

Скорее всего до этого не дойдёт. Где-нибудь при обработке 20 строки числа перестанут помещаться в отведённые для них ячейки. Можно конечно для хранения числителя/знаменателя вместо long int взять long long int = int64_t, тогда это произойдёт попозже, строке так в 40-й.

Цитата Majestio @ 09.06.24, 03:06

amk, ты бы лучше не мне советовал как советовать, а посоветовал бы свой совет ТС'у

Я посоветовал. Воспользоваться не тетрадным вариантом исключения Гаусса, который он реализовал, а нормальным, компьютерным.

[Majestio]

Цитата amk @ 09.06.24, 12:13

Где-нибудь при обработке 20 строки числа перестанут помещаться в отведённые для них ячейки.

Ты не дочитал мой ответ и давай отвечать ... Но ты не на мой текст отвечаешь, а на то, что ты там сам себе придумал! Какие нафик ячейки?!!!

[Qraizer]

H g, тут тема, как я смотрю, сильно переплетена и с инструментами (сиречь Прочим), и с Алгоритмами. Конкретно по Общим вопросам могу лишь предложить double вместо float, а также соптимизировать расчёты, избавивишись от потери точности, где возможно. (Если что, архитектура x86 и некоторых других умеет отлеживать потерю точности при вычисления и в нужным местах кидать эксепшны или хотя бы отражать факт в статусах. Но могут понадобиться инструменты, например, отладчик.) Типа вычитания близких по величине чисел, сложения малых и больших итп. По другим аспектам... в общем, реши, плз, как лучше: отзеркалить, перенести и куда.

Добавлено 09.06.24, 17:35
P.S. Могу порекомендовать boost::multiprecision, когда-то был нужен коллеге, я мне было просто интересно повозиться. За производительность не скажу, но 1000-битные и даже -значные мантиссы легко держала. Интерфейс у неё, правда, не столь удобен, как хотелось бы, так что конкуренты в этом плане у неё есть, думаю.

[Majestio]

Скрытый текст

Цитата Qraizer @ 09.06.24, 17:22

Ребят, ещё немного и в Холивары скину.

Да пусть тут останется. Тема-то не про алгоритм как таковой, а об обеспечение точности.
А холивара не будет. amk, залетел, дан ценные указания использовать GMP и не "велосипедить"
Но произошел обломчег - его либа после двух действий дает 10 == 9 ...

Прикреплённая картинка

Вообще не вижу смысла холиварить

[H g]

Qraizer Если эта тема больше в "Алгоритмы" подходит, можете перенести туда. Правда активность на форуме не столь велика, поэтому перенос в "Алгоритмы" разве что для порядка. По теме : Направления, в которых можно посмотреть, по советам Majestio, amk, Qraizer я получил. Всем спасибо.

Сообщение отредактировано: H g - 10.06.24, 06:12

[H g]

Попробую еще метод Крамера. Определители находятся путем умножения и сложения, а результат, путем деления, но делается это один раз в конце. Результат можно и округлить до нужного знака, в отличие от метота Гаусса где погрешность накапливается от уравнения к уравнению, и "заползает" на результат. Если одно деление, тогда легко можно определить что 9999,99999 это есть 10000.

Так как у меня решение СЛАУ требуется исключительно в прикладных целях, наверное, нужно уточнить пределы задачи. Число уравнений (матрица) максимум до 100 x 100, если будет медленно работать - не критично, такой случай маловероятен. Обычно порядка 20 x 20.
Кроме того матрица сильно разреженная. Вот примерно такая, (решение в MathCAD)

Прикреплённая картинка

Сообщение отредактировано: H g - 10.06.24, 08:40

[MBo]

H g
Итеративные методы пробовали? (например, Гаусса-Зейделя)

[H g]

Цитата MBo @ 10.06.24, 09:00

Итеративные методы пробовали? (например, Гаусса-Зейделя)

Нет, пока с прямыми "воюю" ))
В моем представлении, если ожидается результат какой-то иррациональный, тогда приближаться к нему нужно с помощью итераций.
У меня же наоборот, результат чаще всего целое число, или уж точно из области рациональных чисел.
Я к итерационным способам отношусь с недоверием. Где гарантия того, что на данных, которые введет пользователь метод не улетит черти зная куда, сложнее проверить.

Добавлено 10.06.24, 09:14
Смотрю определители, нужные для метода Крамера хорошо считаются когда много нулей ...
Так может и не плохо получится с моим "велосипедом" )

Сообщение отредактировано: H g - 10.06.24, 09:37

[H g]

Цитата

Попробую еще метод Крамера

Не очень. Если при Гауссе погрешность возникала за счет деления, то в этом методе процесс лезет в область больших чисел. Да и определитель 100 x 100 тоже не подарок. Нужно что то "свое" придумать. Как-то использовать сильную разреженность...

[ЫукпШ]

Цитата H g @ 07.06.24, 21:35

Вопрос как увеличить точность решения ?

В таких случаях я всегда пытаюсь преобразовать
формулу вычислений таким образом, чтобы:
1. деление выполнялось один раз и было последним действием.
даже если это приводит к увеличению количества вычислений
2. произвести последнее деление с "правильным" округлением.
---
Если есть такая возможность, перейти к целым значениям.
Например, если в результате ожидаются дробные миллиметры,
считайте в тысячных долях миллиметров. Главное, чтобы
разрядной сетки целых хватило для всех преобразований.
---
Да, переход к double правильное решение.
Но надо понимать, что использование float/double это
принципиально приближённые вычисления. И если мы начали это
использовать, значит и результат будет приближённый.
Ничего ужасного в этом нет.