Вариационные методы решения СЛАУ (методы минимальных невязок, минимальных поправок, наискорейшего спуска, сопряженных градиентов)

1.

Введение

Постановка задачи.

2. Описание методов решения задачи

2.1 Метод минимальных невязок

2.2 Метод минимальных поправок

2.3 Метод скорейшего спуска

2.4 Метод сопряженных градиентов

3. Алгоритмы и блок-схемы решения

3.1 Метод минимальных невязок

3.2 Метод минимальных поправок

3.3 Метод скорейшего спуска

3.4 Метод сопряженных градиентов

4.Руководство пользователя.

5. Тестирование и оптимизация.

6. Выводы

Список использованных источников

Приложение 1. Контрольные примеры.

Приложение 2. Код программы.

1.

Введение

Постановка задачи

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

Итак, сформулируем задачу:

Дана система линейных уравнений

где A симметричная положительно определенная матрица.

Требуется решить данную СЛАУ следующими итерационными методами вариационного вида:

1. Метод минимальных невязок

2. Метод минимальных поправок

3. Метод скорейшего спуска

4. Метод сопряженных градиентов

Методы решения СЛАУ рассматриваются практически в каждой книге по численным методам. Общее описание алгоритмов для каждого метода дано, например, в книге Самарского А. А., Гулина А. В. «Численные методы». Также популярными книгами являются: М. Ю. Баландин, Э. П. Шурина «Методы решения СЛАУ большой размерности» [1], Y. Saad «Iterative Method for Sparse Linear System». В перечисленных книгах, описаны такие методы как CG (Conjugate Gradient — метод сопряжённых градиентов) и MinRES (Minimal Residual — метод минимальных невязок). В книге изложены не только базовые алгоритмы данных методов, но и указания к их практической реализации.

2. Описание методов решения задачи

2.1 Метод минимальных невязок

Рассмотрим систему с матрицей A=AT>0. Обозначим через

(1)

невязку, которая получается при подстановке приближенного значения xk, полученного на k-й итерации, в уравнение (1). Заметим, что погрешность zk=xkх и невязка rk связаны равенством Azk=rk.

Рассмотрим явный итерационный метод

(2) <=> (3)

где параметр k+l выбирается из условия минимума погрешности ||rk+1|| при заданной норме ||rk||. Метод (2) называется методом минимальных невязок.

Найдем явное выражение для параметра k+l. Из (3) получаем

(4) => (5)

т. е. rk удовлетворяет тому же уравнению, что и погрешность zk=xkх. Возводя обе части уравнения (5) для невязки скалярно в квадрат, получаем

(6)

Отсюда видно, что |||rk+1|| достигает минимума, если

(7)

Таким образом, в методе минимальных невязок переход от k-й итерации к (k+1)-й осуществляется следующим образом. По найденному значению xk вычисляется вектор невязки rk=Axkf и по формуле (7) находится параметр k+l. Затем по формуле (3) досчитывается вектор xk+1.

Метод минимальных невязок (3), (7) сходится с той же скоростью, что и метод простой итерации с оптимальным параметром .

2.2 Метод минимальных поправок

Рассмотрим неявный итерационный метод вида

(8)

Запишем этот метод в виде

. (9)

Вектор k =B1rk называется поправкой на (k+1)-й итерации. Она удовлетворяет тому же уравнению, что и погрешность zk=xkx неявного метода, т. е. уравнению

. (10)

Будем предполагать, что B=BT>0. Методом минимальных поправок называется неявный итерационный метод (8), в котором параметр k+l выбирается из условия минимума нормы ||k+1||B=(Bk+1,k+1)½ при заданном векторе k. В случае B=Е метод минимальных поправок совпадает с методом минимальных невязок.

Найдем выражение для итерационного параметра k+l. Перепишем (10) в виде

и вычислим

.

Отсюда следует, что минимум этого выражения достигается при

. (11)

Для реализации метода минимальных поправок требуется на каждой итерации решить две системы уравнений Bk= rk и Bvk=Ak, откуда найдем вектор vk = B1 Аk (воспользуемся формулой (9)) .

Скорость сходимости метода минимальных поправок определяется границами спектра обобщенной задачи на собственные значения

. (12)

2.3 Метод скорейшего спуска

Рассмотрим явный метод (2) и выберем итерационный параметр k+l из условия минимума ||zk+1||A при заданном векторе zk, где zk+1=xk+1x. Поскольку погрешность zk удовлетворяет уравнению

то

.

Следовательно, это выражение достигает минимума, если

.

Величина zk=xkx здесь неизвестна (так как неизвестно точное решение х), поэтому надо учесть, что Azk=rk=Axkf, и вычислять k+l по формуле

. (13)

2.4 Метод сопряженных градиентов

Пусть A — матрица системы. Рассмотрим следующий класс неявных двухшаговых итерационных методов:

. (14)

Здесь, k+l, k+l — итерационные параметры, которые будут определены далее. Начальное приближение x0 будем задавать произвольно, а вектор x1 вычислять по одношаговой формуле, которая получается из (14) при k=0 и l=1,

. (15)

Если параметры k+l, k+l найдены, то новое приближение xk+l выражается через два предыдущих значения хk и xk1 по формуле

. (16)

Для погрешности zk=xkx из (16) получаем уравнения

.

градиент алгоритм матрица размерность Введем, как и ранее, вспомогательную функцию vk=A½zk, для которой ||vk||=||zk||A. Функция vk удовлетворяет уравнениям

.(17)

где C=A½B1A½. Будем считать, что A и B — симметричные положительно определенные матрицы, удовлетворяющие неравенствам

. (18)

Исключая последовательно векторы vl, v2,, vkl из уравнений (17) находим, что

(19)

где Pk (C) — многочлен степени k от оператора С, удовлетворяющий условию Pk (0)=E.

Поставим задачу выбрать итерационные параметры k, k так, чтобы при любом n=1, 2, … была бы минимальной ||vn||=||zn||A.

Параметр 1 находится из условия минимума ||v1|| при заданном векторе v0. Так же, как и в методе скорейшего спуска, получаем

. (20)

Отметим, что при таком выборе 1 выполняется равенство (Cv1,v0)=0, т. е. векторы v1 и v0 ортогональны в смысле скалярного произведения

.

Далее, рассмотрим погрешность (19) и запишем многочлен Pk© в виде

(21)

где ai (k) — числовые коэффициенты, определяемые параметрами i, i, i=l, 2, …, k. Тогда

. (22)

Найдем условия, которым должны удовлетворять коэффициенты ai (k), минимизирующие ||vn||2. Из (22) получаем

(23)

т. е. ||vn||2 является многочленом второй степени по переменным a1(n), …, an (n). Приравнивая к нулю частные производные ||vn||2/aj (n), j=1, 2, …, n, приходим к системе уравнений

(24)

решение которой ai (n), i=1, 2, …, n и обращает в минимум ||vn||2.

3. Алгоритмы и блок-схемы решения

3.1 Метод минимальных невязок

1. Задаем вектор х0 (начальное приближение).

2. Положим xk = x0, k = 0 (номер итерации)

3. Вычисляем вектор rk = Axk — b (невязка начального приближения).

4. Вычисляем скаляр фk+1 = (rk, Ark) / ||Ark||2.

5. Вычисляем вектор xk+1 = xk + фk+1rk (очередное приближение).

6. Вычисляем вектор rk+1 = rk — фk+1Ark. (невязка k+1 приближения).

7. Проверяем выполнение неравенства ||r (k+1)||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты.

8. Положим k = k + 1 и вернемся к шагу 4.

Рис. 1 Блок-схема метода минимальных невязок

3.2 Метод минимальных поправок

1. Задаем вектор х0 (начальное приближение), матрицу B = BT > 0.

2. Положим xk = x0, k = 0 (номер итерации)

3. Вычисляем вектор rk = Axk — b (невязка начального приближения).

4. Вычисляем вектор щk = rkBk-1.

5. Вычисляем скаляр фk+1 = (Aщk, щk) / (B-1Aщk, Aщk) (шаговый множитель).

6. Вычисляем вектор xk+1 = xk — фk+1 B-1Aщk (очередное приближение).

7. Вычисляем вектор rk+1 = Axk+1 — b. (невязка k+1 приближения).

8. Положим k = k + 1.

9. Проверяем выполнение неравенства ||rk||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты, иначе возвращаемся к шагу 3.

Рис. 2 Блок-схема метода минимальных поправок

3.3 Метод скорейшего спуска

1. Задаем вектор х (0) (начальное приближение).

2. Положим xk = x0, k = 0 (номер итерации)

3. Вычисляем вектор rk = Axk — b (невязка начального приближения).

4. Вычисляем скаляр фk+1 = (rk, rk) / (Ark, Ark) (шаговый множитель).

5. Вычисляем вектор xk+1 = xk — фk+1rk (очередное приближение).

6. Вычисляем вектор rk+1 = Axk+1 — b. (невязка k+1 приближения).

7. Положим k = k + 1.

8. Проверяем выполнение неравенства ||rk||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты, иначе возвращаемся к шагу 3.

Рис. 3 Блок-схема метода скорейшего спуска

3.4 Метод сопряженных градиентов

1. Задаем вектор х (0) (начальное приближение).

2. Положим xk = x0, k = 0 (номер итерации)

3. Вычисляем вектор rk = Axk — b (невязка начального приближения).

4. Положим вектор pk= rk (сопряженный вектор)

5. Вычисляем скаляр бk = (rk, rk) / (Apk, pk).

6. Вычисляем вектор xk+1 = xk + бk pk (очередное приближение)

7. Вычисляем вектор rk+1 = rk — бk Apk. (невязка k+1 приближения)

8. Вычисляем скаляр вk = (rk+1, rk+1) / (rk, rk).

9. Проверяем выполнение неравенства ||rk+1||2 =< 0.001; если «да», то останавливаем работу алгоритма и выводим результаты

10.. Вычисляем pk+1= rk+1+ вkpk

11. Положим k = k + 1.

12. Возвращаемся к шагу 5

Рис. 4 Блок-схема метода сопряжённых градиентов

4.Руководство пользователя

В окне программы пользователю предлагается ввести размерность матрицы системы и требуемую точность, также при желании можно задать величину диагонального преобладания, либо убрать преобладание вообще (см. рис.5)

Рис. 5 Окно программы

При нажатии на кнопку ОК генерируется симметричная матрица заданной размерности, а также правая часть. Далее, при нажатии на кнопку с названием метода, в соответствующий столбец выводится решение и количество итераций (см. рис. 6).

Рис. 6 Решение системы с матрицей размерностью 10×10

Матрицу можно модифицировать, прибавив к главной диагонали заданную величину, для этого нужно ввести необходимую величину в поле «Преобладание», поставить флаг «Модифицировать» и нажать ОК.

Также существует возможность записи сгенерированной матрицы в файл и чтения существующей матрицы из файла. Пользователь может выбрать для записи и чтения один из трёх файлов.

5. Тестирование и оптимизация

Изначально программа тестировалась на симметричных матрицах с диагональным преобладанием, то есть при построении матрицы на главную диагональ прибавлялась 1000, что давало очень хорошую сходимость при заданной точности 0,001 у всех четырёх методов. При этом было видно, что наихудшую скорость сходимости дает метод минимальных поправок, наилучшую — метод сопряжённых градиентов. Разница в скорости сходимости становилась заметна, при размерности матрицы больше чем 200×200 элементов, что видно на графике (см. рис.7)

Рис. 7 Зависимость количества итераций от размерности матрицы.

Для более наглядной демонстрации работы вариационных методов в программу была добавлена возможность генерировать матрицы без диагонального преобладания, а также с диагональным преобладанием, заданным пользователем. В результате были сделаны следующие наблюдения: при отсутствии диагонального преобладания метод минимальных невязок очень медленно сходится (за несколько сотен тысяч итераций), при этом он сходится только для точности не меньше 0,01 и матриц размерностью не больше 50. Аналогичная ситуация наблюдается с методами минимальных поправок и скорейшего спуска. Если задавать небольшое диагональное преобладание, то скорость сходимости увеличивается. На рисунке 8 показано, как изменяется количество итераций в зависимости от диагонального преобладания при размерности матрицы 10×10 и точности, равной 0,01.

Рис. 8 Зависимость количества итераций от диагонального преобладания

Что касается метода сопряжённых градиентов, то диагональное преобладание незначительно влияет на скорость его сходимости, так, в случае матрицы 10×10 элементов без диагонального преобладания и при точности, равной 0,01, метод минимальных невязок сошёлся за 25 329 итераций, а метод сопряжённых градиентов — за 11 итераций.

Рис. 9 Зависимость количества итераций в методе сопряженных градиентов от диагонального преобладания

Таким образом, из всех исследованных методов метод сопряжённых градиентов является наиболее эффективным.

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

6. Выводы

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

Список использованных источников

1. Баландин М. Ю., Шурина Э. П. «Методы решения СЛАУ большой размерности» — Новосибирск: Изд-во НГТУ, 2000

2. Бахвалов Н. С., Жидков Н. М. «Численные методы» — М.: Наука, 1986

3. Самарский А. А., Гулин А. В. «Численные методы», Изд. Наука, 1989

4. Фаддеев В. К., Фаддеева В. Н. «Вычислительные методы линейной алгебры» — М.: Наука 1958.

5. Saad Y. Iterative methods for sparse linear systems. Boston: PWS Publ. Co., 2000.

Приложение 1. Контрольные примеры

Возьмём для наглядности систему с матрицей 5x5 элементов при точности, равной 0,01, и рассмотрим её решения в случаях с различным диагональным преобладанием:

1) Диагональное преобладание отсутствует.

2) На главную диагональ прибавлено 10.

3)На главную диагональ прибавлено 100

4) На главную диагональ прибавлено 1000

Приложение 2. Код программы

//—————————————————————————————————————;

#include

#include

#include

#include

#pragma hdrstop

#include «Unit1.h»

//—————————————————————————————————————;

#pragma package (smart_init)

#pragma resource «*.dfm»

TForm1 *Form1;

int n, m;

float epsilon; //Точность

double**a; //Матрица

double *f; //Столбец правых частей

double *x; //Решение

double *ax;

double**b;

bool ch=0;

double *tm;

//—————————————————————————————————————;

__fastcall TForm1: TForm1(TComponent* Owner)

: TForm (Owner)

{

sg[0]=StringGrid2;

sg[1]=StringGrid3;

sg[2]=StringGrid4;

sg[3]=StringGrid5;

sg[4]=StringGrid6;

sg[5]=StringGrid7;

sg[6]=StringGrid8;

}

//—————————————————————————————————————;

//Матрично-векторное умножение

void mv_mult (double **a, double *x, double *ax)

{

for (int i=0;i

{

ax[i]=0;

for (int j=0;j

{

ax[i]+=a[i][j]*x[j];

}

}

}

//—————————————————————————————————————;

//Векторное вычитание

void vv_substr (double *ax, double *f, double *r)

{

for (int i=0;i

{

r[i]=ax[i]-f[i];

}

}

//—————————————————————————————————————;

//Векторное сложение

void vv_addit (double *ax, double *f, double *r)

{

for (int i=0;i

{

r[i]=ax[i]+f[i];

}

}

//—————————————————————————————————————;

//Норма вектора

double norma (double * vec)

{

double nr=0;

for (int i=0;i

{

nr+=vec[i]*vec[i];

}

return (sqrt (nr));

}

//—————————————————————————————————————;

//Скалярное умножение векторов

double vv_mult (double *a, double *x)

{

double ax=0;

for (int i=0;i

{

ax+=a[i]*x[i];

}

return (ax);

}

//—————————————————————————————————————;

//Проверка на заданную точность нормы

bool check (double *a)

{

double mod;

mod=norma (a);

if (mod

else return (0);

}

//—————————————————————————————————————;

//Построение матрицы B для метода минимальных поправок

void mkB ()

{

ch=1;

b=new double* [n];

for (int i=0; i

{

b[i]=new double[n];

for (int j=0; j

if (i≠j) b[i][j]=0;

else b[i][j]=a[i][j]-300;

}

}

//—————————————————————————————————————;

//Заполнение матрицы

void __fastcall TForm1: Button5Click (TObject *Sender)

{

Button1->Enabled=true;

Button2->Enabled=true;

Button3->Enabled=true;

Button4->Enabled=true;

CheckBox2->Enabled=true;

time_t t;

time (&t); srand ((unsigned)t);

n=StrToInt (Edit1->Text);

m=StrToInt (Edit3->Text);

epsilon=StrToFloat (Edit2->Text);

StringGrid1->ColCount=n;

StringGrid1->RowCount=n;

for (int i=0;i<6;i++)

sg[i]->RowCount=n;

f=new double [n];

if (CheckBox2->Checked==false){

a=new double* [n];

tm=new double [n];

x=new double [n];

for (int i=0;i

a[i]=new double [n];

for (int i=0;i

{

double b=rand ()%200−100;

double c=rand ()%100+1;

x[i]=b/c;

sg[1]->Cells[0][i]=FloatToStr (x[i]);

for (int j=0;j<=i;j++)

{

double b=rand ()%200−100;

double c=rand ()%100+1;

a[i][j]=b/c;

a[j][i]=a[i][j];

if ((i==j)&&(CheckBox1->Checked==true)) a[i][j]+=m;

if (i==j) tm[i]=a[i][j];

StringGrid1->Cells[j][i]=FloatToStr (a[i][j]);

StringGrid1->Cells[i][j]=FloatToStr (a[j][i]);

}

}}

else

for (int i=0;i

{

a[i][i]=tm[i]+(double)m;

tm[i]+=m;

StringGrid1->Cells[i][i]=FloatToStr (a[i][i]);

StringGrid1->Cells[i][i]=FloatToStr (a[i][i]);

}

mv_mult (a, x, f);

for (int i=0;i

sg[0]->Cells[0][i]=FloatToStr (f[i]);

}

//—————————————————————————————————————;

void __fastcall TForm1: FormClose (TObject *Sender, TCloseAction &Action)

{

for (int i=0;i

{

delete [] a[i];

if (ch==1) delete [] b[i];

}

}

//—————————————————————————————————————;

//Метод минимальных невязок

void __fastcall TForm1: Button1Click (TObject *Sender)

{

double *rk; //невязка

double *t; //вспомогательный вектор

double *xk;

double *xk1;

double tau, sm, nrm;

ax=new double [n];

xk=new double [n];

xk1=new double [n];

rk=new double [n];

t=new double [n];

for (int i=0;i

xk[i]=0;

int k=0;

mv_mult (a, xk, ax); //ax=A*xk

vv_substr (ax, f, rk); //Вычисляем невязку: rk=A*xk-f

bool c=check (rk);

while (c==0)

{

mv_mult (a, rk, ax); //A*rk

sm=vv_mult (ax, rk); //Считаем скалярное произведение: (A*rk, rk)

nrm=norma (ax);

tau=sm/(nrm*nrm);

for (int i=0;i

{

t[i]=rk[i]*tau;

ax[i]=ax[i]*tau;

}

vv_substr (xk, t, xk1);

for (int i=0;i

{

xk[i]=xk1[i];

}

vv_substr (rk, ax, t);

for (int i=0;i

{

rk[i]=t[i];

}

c=check (rk);k++;

}

for (int i=0;i

sg[2]->Cells[0][i]=FloatToStr (xk[i]);

sg[6]->Cells[0][0]="Метод мин. невязок" ;

sg[6]->Cells[1][0]=IntToStr (k);;

delete []rk; delete [] xk; delete []xk1; delete [] t;

}

//—————————————————————————————————————;

//Метод минимальных поправок

void __fastcall TForm1: Button2Click (TObject *Sender)

{

double *rk; //невязка

double *wk; //поправка

double *vk;

double *xk;

double *xk1;

double *t;

double tau, sm, smd;

ax=new double [n];

xk=new double [n];

xk1=new double [n];

rk=new double [n];

wk=new double [n];

vk=new double [n];

t=new double [n];

for (int i=0;i

xk[i]=0;

bool c=0;

mkB ();int k=0;

mv_mult (a, xk, ax); //ax=A*xk

vv_substr (ax, f, rk); //Вычисляем невязку: rk=A*xk-f

mv_mult (b, rk, wk); //Вычисляем поправку wk=B*rk

while (c==0)

{

mv_mult (a, wk, ax); //ax=A*wk

mv_mult (b, ax, vk); //vk=B*A*wk

sm=vv_mult (ax, wk); //Считаем скалярное произведение: (A*wk, wk)

smd=vv_mult (vk, ax);

tau=sm/smd;

for (int i=0;i

t[i]=wk[i]*tau;

vv_substr (xk, t, xk1);

for (int i=0;i

{

xk[i]=xk1[i];

vk[i]=vk[i]*tau;

}

vv_substr (wk, vk, t);

for (int i=0;i

{

wk[i]=t[i];

}

c=check (wk);k++;

}

for (int i=0;i

sg[3]->Cells[0][i]=FloatToStr (xk[i]);

sg[6]->Cells[0][1]="Метод мин. поправок" ;

sg[6]->Cells[1][1]=IntToStr (k);

delete []rk; delete [] xk; delete []xk1;delete []wk; delete [] vk; delete [] t;

}

//—————————————————————————————————————;

//Метод скорейшего спуска

void __fastcall TForm1: Button3Click (TObject *Sender)

{

double *rk; //невязка

double *xk;

double *xk1;

double *t;

double tau, sm, smd;

ax=new double [n];

xk=new double [n];

xk1=new double [n];

rk=new double [n];

t=new double [n];

for (int i=0;i

xk[i]=0;

bool c=0;

int k=0;

mv_mult (a, xk, ax); //ax=A*xk

vv_substr (ax, f, rk); //Вычисляем невязку: rk=A*xk-f

while (c==0)

{

mv_mult (a, rk, ax); //A*rk

sm=vv_mult (rk, rk); //Считаем скалярное произведение: (rk, rk)

smd=vv_mult (ax, rk); //Считаем скалярное произведение: (A*rk, rk)

tau=sm/smd;

for (int i=0;i

{

t[i]=rk[i]*tau;

ax[i]=ax[i]*tau;

}

vv_substr (xk, t, xk1);

for (int i=0;i

xk[i]=xk1[i];

vv_substr (rk, ax, t);

for (int i=0;i

{

rk[i]=t[i];

}

c=check (rk);k++;

}

for (int i=0;i

sg[4]->Cells[0][i]=FloatToStr (xk[i]);

sg[6]->Cells[0][2]="Метод наиск. спуска" ;

sg[6]->Cells[1][2]=IntToStr (k);

delete []rk; delete [] xk; delete []xk1;

}

//—————————————————————————————————————;

//Метод сопряженных градиентов

void __fastcall TForm1: Button4Click (TObject *Sender)

{

double *rk; //невязка

double *rk1;

double *pk;

double *xk;

double *xk1;

double *t;

double alpha, beta, sm, smd;

ax=new double [n];

xk=new double [n];

xk1=new double [n];

rk=new double [n];

rk1=new double [n];

pk=new double [n];

t=new double [n];

for (int i=0;i

xk[i]=0;

bool c=0;

int k=0;

mv_mult (a, xk, ax); //ax=A*xk

vv_substr (ax, f, rk); //Вычисляем невязку: r0=A*x0-f

vv_substr (ax, f, pk); //p0=r0

while (c==0)

{

mv_mult (a, pk, ax); //A*pk

sm=vv_mult (rk, rk); //Считаем скалярное произведение: (rk, rk)

smd=vv_mult (ax, pk); //Считаем скалярное произведение: (A*pk, pk)

alpha=sm/smd;

for (int i=0;i

{

t[i]=pk[i]*alpha;

ax[i]=ax[i]*alpha;

}

vv_addit (xk, t, xk1);

for (int i=0;i

xk[i]=xk1[i];

vv_substr (rk, ax, rk1);

sm=vv_mult (rk1,rk1);

smd=vv_mult (rk, rk);

beta=sm/smd;

c=check (rk);

for (int i=0;i

pk[i]=pk[i]*beta;

vv_addit (rk1,pk, pk);

for (int i=0;i

rk[i]=rk1[i];

k++;

}

for (int i=0;i

{

xk[i]=-1*xk[i];

sg[5]->Cells[0][i]=FloatToStr (xk[i]);

}

sg[6]->Cells[0][3]="Метод сопр. градиентов" ;

sg[6]->Cells[1][3]=IntToStr (k);

delete []rk; delete [] xk; delete []xk1;delete []rk1; delete [] pk;

}

//—————————————————————————————————————;

void __fastcall TForm1: Edit1Change (TObject *Sender)

{

CheckBox2->Enabled=false;

CheckBox2->Checked=false;

}

//—————————————————————————————————————;

//Чтение из файла

void __fastcall TForm1: Button6Click (TObject *Sender)

{

FILE *matr;

FILE *pch;

String s=ComboBox1->Text, p=s+" p" ;

Button1->Enabled=true;Button2->Enabled=true;Button3->Enabled=true;

Button4->Enabled=true;CheckBox2->Enabled=true;

epsilon=StrToFloat (Edit2->Text);

matr=fopen (s.c_str ()," rt");

fscanf (matr," %dn" ,&n);

pch=fopen (p.c_str ()," rt");

int tn=StrToInt (Edit1->Text);

if (tn>n) MessageBox (0," Размерность матрицы в файле меньше заданной" ," Ошибка", MB_OK);

else {

n=tn;

f=new double [n];

a=new double* [n];

StringGrid1->ColCount=n;

StringGrid1->RowCount=n;

for (int i=0;i<6;i++)

{

sg[i]->RowCount=n;

sg[1]->Cells[0][i]="" ;

}

for (int i=0;i

a[i]=new double [n];

for (int i=0;i

{

float r;

fscanf (pch," %f" ,&r);

f[i]=r;

sg[0]->Cells[0][i]=FloatToStr (f[i]);

for (int j=0;j

{

fscanf (matr," %f" ,&r);

a[i][j]=r;

StringGrid1->Cells[j][i]=FloatToStr (a[i][j]);

}

}}

fclose (matr);fclose (pch);

}

//—————————————————————————————————————;

//Запись в файл

void __fastcall TForm1: Button7Click (TObject *Sender)

{

FILE *matr;

FILE *pch;

String s=ComboBox2->Text, p=s+" p" ;

matr=fopen (s.c_str ()," w");

fprintf (matr," %dn", n);

pch=fopen (p.c_str ()," w");

for (int i=0;i

{

fprintf (pch," %fn", f[i]);

for (int j=0;j

fprintf (matr," %f «, a[i][j]);

fprintf (matr," n");

}

fclose (matr);fclose (pch);

}