СЛАУ, матрица коэффициентов почти с нулями

[Ben]

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

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

Откуда взялась матрица такого вида? Конечная цель - построить замкнутую поверхность проходящую через заданные точки в трехмерном пространстве. Поверхность задана в сферических координатах, в точках через которые должна проходить поверхность заданы гауссовы функции двух угловых переменных, определяющие радиус вектор $R$ искомой поверхности.
$R=\sum R_1 \exp(-(\theta-\theta_1)^2+(\alpha-\alpha_1)^2)$
Когда углы попадают в заданные значения, гауссиан дает "1", по мере удаления вклад гауссиана быстро убывает, вклад удаленных гауссианов дает значения близкие к нулю. Но это не везде, где-то гауссианы оказываются рядом, там несколько коэффициентов убывает медленно. Чтобы 'притянуть' вклад гауссианов к заданным точкам, решается СЛАУ. Решением СЛАУ являются коэффициенты $R_1$.

[Pphantom]

Ben в сообщении #1440278 писал(а):

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

Есть практически стандартная модификация метода "с выбором ведущего элемента". Идея состоит в том, что при каждом очередном шаге приведения матрицы к треугольному (или диагональному - для модификации Гаусса-Йордана) строки и столбцы переставляются так, чтобы на главной диагонали оказался максимальный по модулю элемент из оставшейся справа и снизу части матрицы. Если переставлять только строки, то это пишется совсем тривиально, если и столбцы - придется учитывать, что при этом порядок элементов в векторе неизвестных меняется (но и с этим можно достаточно легко справиться).

P.S. И не надо набивать формулы долларами внутри. Я поправил ваше сообщение, посмотрите, как это должно выглядеть.

[Ben]

Спасибо! И за ответ, и за доллары!

[arseniiv]

Ну ещё кроме гауссовой элиминации есть и другие методы, связанные часто с разложениями матриц/операторов — правда я не советчик, когда какой уместен, но там не сильно сложно разобраться.

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

[feedinglight]

Если обусловленность очень плохая, то всё плоха.
Если нет, то можно воспользоваться методом Гаусса в форме $LUP$ разложения.
Есть алгоритмы переупорядочивания переменных(Джордж А., Лю Дж. "Численное решение больших разреженных систем уравнений", 1984). Наверно можно просто перенумеровать узлы, чтобы матрица сразу получилась близкой к диагональной (тогда можно использовать просто $LU$ разложение).

Есть итерационные методы с предобусловливанием(BiCGStab, GMRES и т.п.). Раз матрица близка к диагональной, то можно построить предобусловливатель в виде не полного разложения, в профильном формате(строить только значения рядом с диагональю, где значения в исходной матрице достаточно большие) или разреженном, если диагональ не особо выражена. С таким предобусловливателем итерации могут быстро сойтись. Если СЛАУ решается много раз, а матрица не сильно меняется(не совсем понял исходную постановку задачи и метод решения), то может быть будет иметь смысл не на каждой итерации обновлять предобусловливатель. Если элементы исходной матрицы легко вычисляются, то иногда целесообразно вообще не хранить их в памяти.
Это написано просто как бред в свободном стиле.

Как уже отметили, проще всего взять готовый решатель, их куча открытых и бесплатных.

[Евгений Машеров]

Если можно гарантировать, что наибольший коэффициент, равный единице, ровно в одной строке и одном столбце, то проще всего перенумеровать строки и столбцы. А потом решать, как удобно. Гауссом, LU или итеративно.

[ilghiz]

Плюсую на метод с выбором ведущего элемента. Называется по научному LU factorization with partial pivoting, есть готовые алгоритмы в пакете Lapack и работают они очень надежно.

Но если у Вас действительно один элемент в каждой строке существенно больше всех остальных, то подойдет любой итерационный метод (GMRES, BCGStab, QMR, CG c $A^TA$) даже без предобуславливателя. Большинство из этих методов программируется за вечер и код помещается на экран. В википедии есть псевдокод, но если что будет вдруг не понятно, спрашивайте, расскажу.

Достоинства таких методов - при Вашей матрице они сойдутся за очень маленькое (10-20) число итераций на каждой из которых Вам нужно только $N^2$ арифметических операций, в то время, как методы Гаусса - требуют $N^3$ арифметических операций, и от 200 уравнений Вам будет быстрее решать итерационным методом, чем Гауссом.