Кто-нибудь таким маялся? Я, например, реализовал Гауссовский метод, но погрешности умножения и сложения double не дают нормально привести матрицу к треугольной. Вроде бы складываю и умножаю, а вначале 0 - это 0.0000000000000000000001, а через десяток операций это уже 10000000000000. Может кто-нибудь делал/видел алгоритмы. Мне надо решать систему уравнений, например 50х50 как можно быстрее. Причем могут быть нулевые строки. Помогите с алгоритмом...

Как это "нулевые строки"? Такую систему не решить. Может нулевые диагональные элементы? Тогда придеться делать соответствующую проверку и менять строки местами. А вообще лучше использовать метод Ньютона или Зейделя. Имхо быстрее и проще.

А исходники есть? Нулевые строки, я имел в виду, что есть 0 x1+0 x2+0 x3+0 x4=0. Я их могу вычеркнуть, но тогда матрица не будет квадратная

Помню, я, учась в университете, реализовывал пару-тройку таких алгоритмов.
Даже сравнение количества операций было.
Вот тока не помню, остались ли у меня такие доки. По-моему, сами реализации где-то валялись.
Приду домой - гляну.
Тебе обязательно на С++? По-моему, там у меня Borland C++ && Delphi.
Накрайняк, если ничего хорошего не найдешь и я у себя в эл. виде не найду - напишу что-то типа "статьи" - хотя, вроде, в и-нете полно таких доков..
Кстати.. ИМХО, любой алгоритм будет работать с погрешностью - ведь комп обрабатывает числа с плавающей точкой с некоторой степенью точности.
Я реализовывал алгоритмы, к-е теоретически должны давать точный результат, но на практике появлялись погрешности.
Тебе что важнее? Скорость или точность?
Хорошую скорость дают методы Зейделя, метод какого-то-там ( названия не помню, но маялся я с ним долго..) спуска... короче, надо добраться домой.

Если проблема - только с точностью, то используй библиотеки высокой точности для промеж. вычислений:
1) набирай коэфф. уравнений в форме double.
2) для вычислений переведи их в какой-нибудь более точный вид данных (в соответствии с библиотекой высокой точности),
3) считай что тебе нужно, а
4) результат обратно переводи ("downcast") в double

теперь вопрос какой метод решения выбирать и в какой форме хранить:
если твои матрицы будут очень жидконаселёнными, стоит хранить в форме
массив списков а не 2-х мерный массив.
а метод решения ты уж сам выбирай, хоть Kramer , хоть Gauss, хоть Seidel-Gauss... как тебе угодно.

Цитата

А исходники есть? Нулевые строки, я имел в виду, что есть 0 x1+0 x2+0 x3+0 x4=0. Я их могу вычеркнуть, но тогда матрица не будет квадратная

Эй-эй-эй!
а ты уверен, что твои системы разрешимы и меют тока ОДНО решение?
По-моему, все вышеназванные методы, тока для системы, у к-й одно и только одно решение!

Если проблема - только с точностью, то используй библиотеки высокой точности для промеж. вычислений:
1) набирай коэфф. уравнений в форме double.
2) для вычислений переведи их в какой-нибудь более точный вид данных (в соответствии с библиотекой высокой точности),
3) считай что тебе нужно, а
4) результат обратно переводи ("downcast") в double

теперь вопрос какой метод решения выбирать и в какой форме хранить:
если твои матрицы будут очень жидконаселёнными, стоит хранить в форме
массив списков а не 2-х мерный массив.
а метод решения ты уж сам выбирай, хоть Kramer , хоть Gauss, хоть Seidel-Gauss... как тебе угодно.

Вот алгоритм приведения матрицы к диагональному виду(написано для вычисления определителя):

//Приведение матрицы к треугольному виду
int fn(double **a, int size)
{
for(int k=0; k<size-1; k++)//Прогон по диагонали
{
for(int j=k+1;j<size;j++)
{
double koef=(a[j][k]/a[k][k]);
a[j][k]=0;
for(int i=k+1;i<size;i++)
{
a[j][i]+=(-koef)*a[k][i];
}
}
}
}
Если вдруг будут нулевые диагональные элементы (a[k][k]=0), то алгоритм не прокатит, нужно ввести дополнительное условие типа:
if(a[k][k]==0){меняем местами текущую и следующую строки}

Проверял результат (определитель матрицы) в Матлабе - результаты идентичны до восьмого знака после запятой (дальше не смотрел). Так что дабл вроде канает. Считает правда дольше чем Матлаб. Матрицу 150х150 на 1500ГГц считает где-то пол минуты.

Если нужен "умный" решатель, надо заранее просмотреть систему на "квадратность", нулевые строки, сингулярность, а уж потом решать любым из методов.

Мне советовали использовать метод Холецкого. Но у меня не симметричная матрица. Как он будет работать для таковой?

Слушай, у меня есть моя отличная реализация этой задачи на Delphi.
Хочешь я тебе её подгоню, а ты там сам в коде поковыряешься - там отличные алгоритмы.
Ну так как?

Цитата

int gauss(double *aMatr,double *bMatr,double *vector,int razmer); void copy_to_temp(double *pVar,double *pTemp,int n); void umn_temp(double *pTemp,double x,int n); void plus_temp(double *pVar,double *pTemp,int n); int poisk_ne_null(int m,double *pVar,int n); void zamena(double *pVar1,double *pVar2,int n); void zamena(double *pVar1,double *pVar2,int n) { double t; n++; for(int i=0;i<n;i++) {   t=*(pVar1+i);   *(pVar1+i)=*(pVar2+i);   *(pVar2+i)=t; } } int poisk_ne_null(int m,double *pVar,int n) { int k=0,i; for(i=m;i<n;i++) {   if(*(pVar+k*100))   return i;   k++; } return -1; } void plus_temp(double *pVar,double *pTemp,int n) { n++; for(int i=0;i<n;i++) {   *(pVar+i)+=*(pTemp+i); } } void umn_temp(double *pTemp,double x,int n) { n++; for(int i=0;i<n;i++) {   *(pTemp+i)=x*(*(pTemp+i)); } } int gauss(double *aMatr,double *bMatr,double *vector,int razmer) { double A[20][20],XX[20],temp[20]; int i,j; //--------------------------------------------------------- for(i=0;i<razmer;i++) {   for(j=0;j<razmer;j++)   {   A[i][j]=*(aMatr+i*razmer+j);   } } for(i=0;i<razmer;i++) {   A[i][razmer]=*(bMatr+i); } //----------------------------------------------------------- for(i=0;i<razmer-1;i++) {   if(A[i][i]==0)   {   j=poisk_ne_null(i,&A[i][i],razmer);   if(j==-1)     return 0;   else     zamena(&A[i][0],&A[j][0],razmer);   }   for(j=i+1;j<razmer;j++)   {   copy_to_temp(&A[i][0],&temp[0],razmer);   umn_temp(&temp[0],-A[j][i]/A[i][i],razmer);   plus_temp(&A[j][0],&temp[0],razmer);   } } for(i=razmer-1;i>0;i--) {   if(A[i][i]==0)   return 0;   for(j=i-1;j>=0;j--)   {   copy_to_temp(&A[i][0],&temp[0],razmer);   umn_temp(&temp[0],-A[j][i]/A[i][i],razmer);   plus_temp(&A[j][0],&temp[0],razmer);   } } for(i=0;i<razmer;i++) {   XX[i]=A[i][razmer]/A[i][i]; } //----------------------------------------------------------- for(i=0;i<razmer;i++) {   *(vector+i)=XX[i]; } //----------------------------------------------------------- return 1; } void copy_to_temp(double *pVar,double *pTemp,int n) { n++; for(int i=0;i<n;i++) {   *(pTemp+i)=*(pVar+i); } }

Можно многое оптимизировать, но это нормальная рабочая версия...

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

Если решать систему с линейно зависимыми строками то в итоге придется вырпжать все через некоторые переменные

Пример:

Цитата

X1+X2=5; 2*X1+2*X2=10;

Тогда решение выглядит так:

Цитата

X1=5-X2;

Те все это придется выражать через X2...

Так у тебя что задача написать прогу которая решает любые системы??? Если так то давай займемся...Без б, помогу...

Решение систем линейных алгебраических уравнений (СЛАУ) в общем виде непростая задача.
Похоже у Вас плохо обусловленная или вырожденная СЛАУ. Для ее решения надо применять специальные устойчивые методы регуляризации или сингулярного разложения.
Ощий подход сдедующий:
1) решаем методом Гаусса (он самый быстрый), если не решает, то
2) решаем либо методом сингулярного разложения (svd) (сложный, медленный, но надежный),
3) либо методом наискорейшего градиентного спуска (простой, надежный, но медленный).
Для более точных рекомендаций надо уточнить цель Ваших расчетов.
Если проблема еще не решена, то пишите, я вышлю нужные процедуры.
Н.И. Клевец
[email protected]

Ребят, я тут новенький!))) может у кого есть готовый алгоритм решения СЛАУ методом наискорейшего градиентного спуска?! Если есть и нетрудно - вышлите пожалуйста - он мне ну очень срочно нужен. Заранее спасибо!

[email protected]