Форум программистов, компьютерный форум, киберфорум
С++ для начинающих
Войти
Регистрация
Восстановить пароль
 
Рейтинг 4.81/27: Рейтинг темы: голосов - 27, средняя оценка - 4.81
46 / 46 / 14
Регистрация: 30.03.2009
Сообщений: 523
1

Решение СЛАУ большой размерности методом сопряженных градиентов

19.12.2012, 17:28. Просмотров 5095. Ответов 3

Всем првиет! Возникла проблемка с методом сопряженных градиентов. Если задавать самому значения матрицы и правой части, то все решается отлично. Сделал rand'омное заполнение. В итоге получаю, что программа работает только с матрицами размера 2*2 и то не со всеми, даже 3*3 она уже не считает правильно, в результате выводит ответ в виде -1.#IND. А мне нужно, чтобы она работала с матрицами большого порядка, мне надо узнать время выполнения алгоритма. В чем можнт быть проблема? Я даже double заменил на long double и все равно не помогает.
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
#include<stdio.h> 
#include <time.h>
#include <conio.h>
#include <iostream>
using namespace std;
 
#define M 2 //размерность системы
 
#define E 0.1 //точность вычислений
 
 
long double Xk[M], Zk[M];
long double Rk[M], Sz[M], alpha, beta, mf;
long double Spr, Spr1, Spz;
long double A[M][M];//={2.0, -1.0, 0.0, -1.0, 5.0, 3.0, 0.0, 3.0, 1.0};
long double F[M] ;//= {1.0, 7.0, 4.0};
 
int main(void){
    int  i, j, kl = 1;
 
    srand(time(0));
    for (i = 0; i < M; i++)
    {
        for (j = 0; j < M; j++){
     A[i][j] = rand() % 2 - rand() % 1; // заполняем массив случайными значениями в диапазоне от -1 до 1 включительно
     cout << A[i][j] << " "; // печать элементов одномерного массива
        }
        cout << endl;
    }
 
    for (j = 0; j < M; j++){
     F[j] = rand() % 2 - rand() % 2; // заполняем массив случайными значениями в диапазоне от -1 до 1 включительно
     cout << F[j] << " "; // печать элементов одномерного массива 
        }
 
    clock_t time;
    time = clock();//начальное время
 
/* Вычисляем сумму квадратов элементов вектора F*/  
    for(mf = 0,i = 0; i < M; i++){
        mf += F[i] * F[i];
      }
 
 
/* Задаем начальное приближение корней. В Хk хранятся значения корней
 * к-й итерации. */
    for(i = 0; i < M; i++){
        Xk[i] = 0.2;
    } 
 
/* Задаем начальное значение r0 и z0. */
    for(i = 0; i < M; i++){
        for(Sz[i]=0,j = 0; j < M; j++)
            Sz[i] += A[i][j] * Xk[j];
      Rk[i] = F[i] - Sz[i];
      Zk[i] = Rk[i];  
      }
    
    do{ 
  /* Вычисляем числитель и знаменатель для коэффициента
   * alpha = (rk-1,rk-1)/(Azk-1,zk-1) */
        Spz = 0;
        Spr = 0;
        for(i = 0; i < M; i++){ 
            for(Sz[i] = 0, j = 0; j < M; j++){
                Sz[i] += A[i][j] * Zk[j];
            }
          Spz += Sz[i] * Zk[i];
          Spr += Rk[i] * Rk[i];
          }
        alpha = Spr / Spz;             /*  alpha    */
    
 
/* Вычисляем вектор решения: xk = xk-1+ alpha * zk-1, 
    вектор невязки: rk = rk-1 - alpha * A * zk-1 и числитель для betaa равный (rk,rk) */
        Spr1 = 0;
        for(i = 0; i < M; i++){
            Xk[i] += alpha * Zk[i];
            Rk[i] -= alpha * Sz[i];
            Spr1 += Rk[i] * Rk[i];
            cout << "Iter #" << kl;
            cout << " " << "X[" << i << "] = " << Xk[i] << endl;
        }
        cout << endl;
        kl++;
 
/* Вычисляем  beta  */
        beta = Spr1 / Spr; 
 
/* Вычисляем вектор спуска: zk = rk+ beta * zk-1 */
        for(i = 0; i < M; i++)
            Zk[i] = Rk[i] + beta * Zk[i];       
      }
/* Проверяем условие выхода из итерационного цикла  */
    while(Spr1 / mf > E * E);
 
    time = clock() - time;//время выполнения
    cout << "Time result = " << (double)time << "ms" << endl;
    system("pause");
    return 0;
  }
Добавлено через 59 минут
UPD.Понял почему не считает! Матрица же симметричная и положительно определенная должна быть!!! Теперь возникает другой вопрос: как рандомно заполнить симметричную матрицу?

Добавлено через 1 час 1 минуту
UPD.Все. Во всем разобрался
0
Programming
Эксперт
94731 / 64177 / 26122
Регистрация: 12.04.2006
Сообщений: 116,782
19.12.2012, 17:28
Ответы с готовыми решениями:

Решение СЛАУ методом сопряженных градиентов (нужна литература)
Подскажите пожалуйста литературу, в которой подробно описывается метод решения СЛАУ методом...

СЛАУ методом сопряженных градиентов Delphi
Всем привет. Хочу заказать работу. Кто может помочь? Литература, алгоритм есть.

Метод сопряженных градиентов для решения СЛАУ
Пишу прогу для решения СЛАУ методом сопряженных градиентов. Есть СЛАУ, Ax=B, заданная матрицами ...

Минимизация методом сопряженных градиентов
Помогите, пожалуйста, реализовать нахождение минимума функции методом сопряженных градиентов в...

3
0 / 0 / 0
Регистрация: 25.12.2010
Сообщений: 4
02.05.2013, 17:39 2
}{@k, не могли бы вы скинуть своё решение, пожалуйста) очень нужно
0
46 / 46 / 14
Регистрация: 30.03.2009
Сообщений: 523
05.05.2013, 21:15  [ТС] 3
аия, держите, только там еще и метод Гаусса реализован для сравнения скорости вычисления, но его несложно удалить
C++
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
#include <stdio.h> 
#include <ctime>
#include <conio.h>
#include <iostream>
using namespace std;
 
#define M 50 //размерность системы
#define E 0.01 //точность вычислений
 
float Xk[M], Zk[M];
float Rk[M], Sz[M], alpha, beta, mf;
float Spr, Spr1, Spz;
float A[M][M];
float F[M] ;
 
float conj_grad(float A[M][M], float F[M]){
    int i,j,kl=1;
    double max_iter = 100000;
/* Вычисляем сумму квадратов элементов вектора F*/  
    for(mf = 0,i = 0; i < M; i++){
        mf += F[i] * F[i];
      }
 
 
/* Задаем начальное приближение корней. В Хk хранятся значения корней
 * к-й итерации. */
    for(i = 0; i < M; i++){
        Xk[i] = 0.2;
    } 
 
/* Задаем начальное значение r0 и z0. */
    for(i = 0; i < M; i++){
        for(Sz[i]=0,j = 0; j < M; j++)
            Sz[i] += A[i][j] * Xk[j];
      Rk[i] = F[i] - Sz[i];
      Zk[i] = Rk[i];  
      }
    int Iteration = 0;
    do{ 
        Iteration++;
  /* Вычисляем числитель и знаменатель для коэффициента
   * alpha = (rk-1,rk-1)/(Azk-1,zk-1) */
        Spz = 0;
        Spr = 0;
        for(i = 0; i < M; i++){ 
            for(Sz[i] = 0, j = 0; j < M; j++){
                Sz[i] += A[i][j] * Zk[j];
            }
          Spz += Sz[i] * Zk[i];
          Spr += Rk[i] * Rk[i];
          }
        alpha = Spr / Spz;             /*  alpha    */
    
 
/* Вычисляем вектор решения: xk = xk-1+ alpha * zk-1, 
    вектор невязки: rk = rk-1 - alpha * A * zk-1 и числитель для betaa равный (rk,rk) */
        Spr1 = 0;
        for(i = 0; i < M; i++){
            Xk[i] += alpha * Zk[i];
            Rk[i] -= alpha * Sz[i];
            Spr1 += Rk[i] * Rk[i];
            //cout << "Iter #" << kl;
            //cout << " " << "X[" << i << "] = " << Xk[i] << endl;
        }
        //cout << endl;
        kl++;
        cout << "Rko-vo iter" << kl << endl;
 
/* Вычисляем  beta  */
        beta = Spr1 / Spr; 
 
/* Вычисляем вектор спуска: zk = rk+ beta * zk-1 */
        for(i = 0; i < M; i++)
            Zk[i] = Rk[i] + beta * Zk[i];       
      }
/* Проверяем условие выхода из итерационного цикла  */
    while(Spr1 / mf > E * E && Iteration < max_iter);
 
    return 0;
}
 
void GlavElem(int k, float mas[][M + 1], int otv[])
{
  int i, j, i_max = k, j_max = k;
  float temp;
  //Ищем максимальный по модулю элемент
  for (i = k; i < M; i++)
  {
      for (j = k; j < M; j++)
      {
          if (fabs(mas[i_max][j_max]) < fabs(mas[i][j]))
          {
              i_max = i;
              j_max = j;
          }
      }
  }
  //Переставляем строки
  for (j = k; j < M + 1; j++)
  {
      temp = mas[k][j];
      mas[k][j] = mas[i_max][j];
      mas[i_max][j] = temp;
  }
  //Переставляем столбцы
  for (i = 0; i < M; i++)
  {
      temp = mas[i] [k];
      mas[i] [k] = mas[i] [j_max];
      mas[i] [j_max] = temp;
  }
  //Учитываем изменение порядка корней
  i = otv[k];
  otv[k] = otv[j_max];
  otv[j_max] = i;
}
 
float gauss(float A[M][M], float F[M]){
    int i, j, k; int otv[M];
    float x[M], mas[M][M + 1];
 
    for (i = 0; i < M; i++){
        for (j = 0; j < M; j++){
        mas[i][j] = A[i][j];
        }
        mas[i][M] = F[i];
    }
    
    for (i = 0; i < M; i++)
  {
      otv[i] = i;
  }
    for (k = 0; k < M; k++)
  { 
      //На какой позиции должен стоять главный элемент
      GlavElem(k, mas, otv); //Установка главного элемента
          
 
    for (j = M; j >= k; j--)
    {
        mas[k][j] /= mas[k][k];
    }
 
    for (i = k + 1; i < M; i++)
    {
        for (j = M; j >= k; j--)
        {
            mas[i][j] -= mas[k][j] * mas[i][k];
        }
    }
  }
  //Обратный ход
  for (i = 0; i < M; i++)
  {
      x[i] = mas[i][M];
  }
  for (i = M - 2; i >= 0; i--)
  {
      for (j = i + 1; j < M; j++)
      {
          x[i] -= x[j] * mas[i][j];
      }
  }
 
    return 0;
}
 
int main(void){
    int  i, j, kl = 1;
 
    srand(time(0));
    
    for (i = 0; i < M; i++){
        for (j = 0; j < M; j++){
            A[i][i] = rand() % 2;
            if ((i != j) && (i < j)){
                A[i][j] = rand() % 2 - rand() % 2;//заполняем массив случайными значениями
            }
                A[j][i] = A[i][j];
        }
    }
 
    for (j = 0; j < M; j++){
     F[j] = rand() % 2 - rand() % 2; // заполняем массив случайными значениями
    }
 
    clock_t time_conj_grad; 
    clock_t time_gauss;
    
    time_conj_grad = clock();//начальное время    
    conj_grad(A,F); 
    time_conj_grad = clock() - time_conj_grad;//время выполнения
    cout << "Running time the conjugate gradient algorithm = " << (double)time_conj_grad << "ms" << endl;
    
    cout << endl <<endl;
 
    time_gauss = clock();//начальное время
    gauss(A,F);
    time_gauss = clock() - time_gauss;//время выполнения
    cout << "Running time Gauss = " << (double)time_gauss << "ms" << endl << endl;
    system("pause");
    return 0;
  }
0
0 / 0 / 0
Регистрация: 25.12.2010
Сообщений: 4
06.05.2013, 13:11 4
}{@k, спасибо!
0
IT_Exp
Эксперт
87844 / 49110 / 22898
Регистрация: 17.06.2006
Сообщений: 92,604
06.05.2013, 13:11

Заказываю контрольные, курсовые, дипломные и любые другие студенческие работы здесь.

Метод сопряжённых градиентов
Всем доброго времени суток. Хочу разобраться в методе сопряжённых градиентов. Нашёл в сети задачу...

Метод сопряженных градиентов
Где в excel 2010 метод сопряженных градиентов (поиск решения)?

Метод Сопряжённых Градиентов
#include &lt;iostream&gt; #include &lt;cstdlib&gt; #include &lt;cmath&gt; #include &quot;sol.h&quot; using namespace std;...

метод сопряжённых градиентов
Решение задачи оптимизации с помощью метода сопряжённых градиентов существенно зависит от...


Искать еще темы с ответами

Или воспользуйтесь поиском по форуму:
4
Ответ Создать тему
Опции темы

КиберФорум - форум программистов, компьютерный форум, программирование
Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.