99703

Программа реализующяа решение плохо обусловленных систем методом регуляции А.Н. Тихонова

Курсовая

Информатика, кибернетика и программирование

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

Русский

2016-10-08

117.26 KB

6 чел.

Содержание

  1.  Постановка задачи………………………………………………………….…..1
  2.  Теоретические сведения……………………………………………………….2
  3.  Реализация алгоритма Тихонова……………………………………………...8
  4.  Руководство пользователя……………………………………………………15
  5.  Вывод…………………………………………………………………………..15

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

Написать программу реализующую решение плохо обусловленных систем методом регуляции А.Н.Тихонова

Теоретические сведения.

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

Пусть имеем систему линейных алгебраических уравнений  Ax=b (1.1). Предположим, что матрице А задана, является не вырожденной и плохо обусловленной.  В  этом  случае  решение  системы (1.1) существует, единственно  и  неустойчиво,  т.е.  задача  некорректна.  Некорректно поставленные  задачи  создают  значительные  трудности  при  их приближенном решении. Как правило, исходные данные прикладных задач известны с некоторой погрешностью. Отыскание неустойчивых решений прикладных  задач  приводит  к  тому,  что  классические  методы вычислительной  математики  могут дать  результат,  сколь  угодно  сильно отличающийся  от точного решения,  при условии, что  исходные данные отличаются сколь угодно мало.

Рассмотри простейший случай, когда матрица А - симметричная матрица. Пусть  –  ее собственные значения, упорядоченные в порядке убывания . Соответствующую ортогональную систему собственных векторов обозначим через   . Решением  системы Аx=b является вектор

X=, где .

При реально заданной правой части  решением будет

.

Коэффициент  может оказаться очень большим при малых ,  что ведет к сильному искажению решения. Иногда заранее известно, что в разложении искомого решения задачи  коэффициенты сi  , соответствующие малым по модулю , малы. В этом случае следует принять какие-то  меры с тем, что бы «отфильтровать»  эти составляющие решения.

Рассмотрим способы, иллюстрированные на примерах матриц с >0.

Первый способ.

Задавая некоторое α>0, находим решение xα системы

(A+ αΕ)xα =.

Оно записывается в виде xα = . Так как   , то наличие малого параметра  не существенно изменит слагаемое, соответствующие  большим . В тоже время при  имеем <<. это означает что введение параметра   приводит к существенному уменьшению роли слагаемых, соответствующим малым . Подбор оптимального значения  обычно осуществляется экспериментально, сравнивая результаты расчетов при различных .

Второй способ.

Заключается в следующем. Будем решать систему уравнений каким-то итерационным способом. Рассмотрим случай итерации по формуле

(*)

При некотором начальном приближении x0. Пусть

= ,  

Подставим эти выражения в (*) и, приравнивая коэффициенты  про ei , получим соотношение .

Последовательно подставляя  через предыдущие, имеем  =.

Если <1 , то при n 

 

Поэтому . Пусть  При больших значениях      быстро стремится к 0 с ростом n и   близко к своему предельному значению.

В других случаях решение задачи находят, минимизируя некоторый функционал, близкий к F(x)=(Ax,x)-2(b,x), например функционал F(x)+ (x,x)  с малым

Успешность  применения описанных приемов в случае несимметричных матриц А в существенной степени  зависит от  структуры жордановой формы и от ряда свойств матрицы. Здесь часто решение находят минимизируя функционал (Аx-b, Аx-b)+(x,x), при малых >0.

Метод А.Н.Тихонова применим к задаче (1.1)  в ее общей постановке (без условия самосопряженности).   

В  общей  постановке  некорректной  задачи (1.1)  рассмотрим  как совместную  неопределенную  систему  линейных  уравнений. Такая система имеет бесконечно много решений, поэтому надо ввести дополнительные условия, чтобы задача стала определенной.  

Во многих прикладных задачах таким условием является требование, чтобы длина вектора решения х была минимальной среди всех решений этой системы, т.е. чтобы скалярное произведение (х,х) было минимальным. Такое решение системы (1.1) называется нормальным и обозначается xn . Вычисление на ЭВМ минора матрицы А, равного нулю, обычно приводит к результату отличному, от нуля даже при малых ошибках округления. В такой ситуации совершенно не ясно, какой минор нужно считать базисным минором  матрицы А. Алгоритм становится некорректным.

Для  решения  некорректных  математических  задач  А.Н.Тихоновым  был предложен  метод  регуляризации.  Применительно  к  нашей  задаче  метод регуляризации основан на решении не исходной задачи, которая является неустойчивой (или  неопределенной),  а  близкой  к  ней  вспомогательной, содержащей некоторый положительный параметр регуляризации α.  После  этого  новая  задача  оказывается  устойчивой (или  однозначно разрешимой)  для  любого  выбора  положительного  параметра  α . Предположим,  что  в  линейном  пространстве  Rn  введено  скалярное произведение, перепишем систему (1.1) в эквивалентной форме:

x-b, Аx-b) = 0. (1.2)  

Регуляризованная  задача  формулируется  так:

Найти  решение  xα, зависящее от α , при котором функционал Fα(x) =(Ax-b, Ax-b)+α(x,x)  (1.3), где  α>0 , является  минимальным,  и  выбрать  параметр  α    так,  чтобы регуляризованное  значение хα  было  наиболее  близким  к  нормальному решению хn системы (1.1) . Регуляризирующий функционал Fα(x) запишем, вычисляя скалярное произведение в ортонормированном базисе:  

Fα(x) =     (1.4)  

Должно выполняться необходимое условие минимума функции:

Fα(x)=Fα(…, , поэтому имеем

=, j=   (1.5)

Равенство (1.5) в матричной форме имеет вид:

(ATA+ αΕ)xα =AT b,α>0     (1.6)

Таким  образом,  решение  неопределенной  системы (1.1) сведено  к решению однозначно разрешимой системы (1.6). Решение этой  системы можно  найти,  используя  прямые  или  итерационные  метода  решения корректных задач. Таким образом, простейший вариант метода А.Н. Тихонова заключается в отыскании точки минимума функционала  

Fα(x)= α

x      (1.7)

где α - малый положительный параметр.

Эта задача на минимум равносильно определению x из уравнения   

?????????αx+Α*Ax=A*b.

Пример.

Для реализации численного примера был выбран метод Тихонова решения плохо обусловленных СЛАУ.  В качестве исходной была взята СЛАУ Ax=b, имеющая в матричной записи вид:

Определитель матрицы коэффициентов этой системы близок к нулю – он равен 0.000125. Попробуем решить эту систему с помощью обратной матрицы:

x=A-1b

Получим    x1=316

                  x2=-990

                  x3=832

Теперь предположим, что правая часть нам известна приближенно, с погрешностью 0.1  Изменим, к примеру, третий элемент вектора-столбца с 1 на 1.1 :

Попробуем решить новую систему также с помощью обратного оператора. Мы получаем другой результат:

               x1=348

       x2=-1090

       x3=916.

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

Будем искать решение методом Тихонова. В теоретической части было показано, что целесообразно использовать регуляризирующий оператор следующего вида:   (ATA+ αΕ)xα =AT b, где E – единичная матрица, xα  – приближенное нормальное решение, AT – транспонированная исходная матрица, α – параметр регуляризации,       b – правая часть, заданная неточно. Эту задачу можно решать стандартными методами, задав предварительно функцию α = α (b), удовлетворяющую условиям теоремы Тихонова.

Приближение к нормальному решению

X(1)= 3.47

X(2)=-1.08

X(3)= 9.15

Значение правой части при подстановке прибл. решения

B1(1)= 9.99

B1(2)= 1.00

B1(3)= 1.09

Реализация алгоритма Тихонова.

Создаем транспонированную матрицу _At матрицы _A[j][i], которая, в свою очередь является главной матрицей программы

void _TRANS_A(){                 

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

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

    _At[i][j]= _A[j][i];

}

//-----------------------------------------------------

Создаем единичную матрицу   _ONE[i][j]

for(int i=0;i<size - 1;i++)

 BB[i] = new double [size];

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

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

  if(i==j)

    _ONE[i][j] = 1;

   else

                _ONE[i][j] = 0;

//-----------------------------------------------------

Определяем знак коэффициента при малом определителе

double _COEFF(int _FPARAM,int _SPARAM){

  int _SUMM = _FPARAM + _SPARAM;

  double _TMP = 1.0;

     for(int k=0;k<_SUMM;k++)

   _TMP*=(-1);

  return _TMP;

}

//-----------------------------------------------------

считает определитель на 1 размер меньше исходной матрицы

double _ALL_ADD(int _FPARAM,int _SPARAM){

         double **TMP = new double*[size - 1];

    for(int i=0;i<size-1;i++)

   TMP[i] = new double [size - 1];

    for(int i=0;i<size - 1;i++)

     for(int j=0;j<size - 1;j++)

      TMP[i][j] = 0;

       for(int v=0;v<size;v++)

     for(int w=0;w<size;w++)

      {

        AA[v][w]._VALUE = _A2[v][w];

        AA[v][w]._FLAG = true;

      }

     for(int v=0;v<size;v++)

      {

              AA[_FPARAM][v]._FLAG = false;//фолсим строку 

              AA[v][_SPARAM]._FLAG = false;//фолсим столбец

      }

    k = 0;

    int x = 0,y=0;

     {

    for(int v=0;v<size;v++)

     for(int w=0;w<size;w++)

      {

      if(AA[v][w]._FLAG)

       {

          TMP[x][y] = AA[v][w]._VALUE;

          y++;

           if(y == size - 1)

         {

           x++;

           y = 0;

         }

        if(x == size - 1) {break;break;}

       }

      }

     }

   

    double res = DET(TMP,size - 1);

    delete [] TMP;

    return res;

    

}

//------------------------------------------------------------

определитель матрицы А2; где А2==АТ *А+α*Е (для разложения используем первую строчку)

double DET_COUNT(){

 double S1;

  for(int z=0;z<size;z++)

   S1+=_A2[0][z]*_COEFF(0,z)*_ALL_ADD(0,z);

 return S1;

 }

//-------------------------------------------------------------

Находим обратную матрицу А2==_Ar[j][i]

void _REV_MATR(){

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

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

  _Ar[j][i] = (double)(_COEFF(i,j)*_ALL_ADD(i,j)/DET_COUNT());

 }

//-------------------------------------------------------------

Определяем разницу между правой и левой частью уравнения по абсолютному значению и сравниваем с ε, т.е. сравниваем |АТ*В- (АТ *А+α*Е)*Хα|<ε, если условие выполниется, то заканчиваем цикл

bool _ALL_RIGHT(){   

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

 {

    if(abs(RIGHT[i] - _U1[i]) > precision)

    return false;

 }

 return true;

}

//-------------------------------------------------------------

Функция, результатом которой является _FPARAM в степени _SPARAM

double _POW(double _FPARAM,int _SPARAM){

  double S2 = 1;

  int z;

       if(_SPARAM==0)

  return 1;

 else

        for(z=0;z<_SPARAM;z++)

   S2*=_FPARAM;

   return S2;

}

//------------------------------------------------------------

void MAIN(){

 обозначаем, что _A[i][j] главная  матрица

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

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

_A[i][j] = MAIN_MATRIX[i][j];

транспонируем главную матрицу

   _TRANS_A();      

                                            

транспонированная матрицу _At[i][j] умножаем на правый столбец (получаем правую часть уравнения Тихонова==_B[i])

   for(int i=0;i<size;i++){      

     _S = 0;                    

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

    {

        _S+=_At[i][j]*RIGHT[j];

        _B[i] = _S;

     

       }

    }

Транспонированную главную матрицу умножаем на главную матрицу получем матрицу_A1[i][j], т.е. _A1[i][j]==АТ

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

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

    {

    _S = 0;

    for(int k=0;k<size;k++){

     _S+= _At[i][k]*_A[k][j];

         _A1[i][j] = _S;

     }

    }

     Сначала количесвто итераций равно 0

             q = 0;

    do{

рассчитываем α как количество итераций деленное 4 в  степени количества итерций, т.е находим параметр регуляризации α

     this->_ALPHA = (double)(q / _POW(4,q));

    

находим матрицу равную произведению α и единичной матрицы, т.е. _One1[i][j]==α*Е

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

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

      _One1[i][j] = _ONE[i][j]*_ALPHA;

находим матрицу  _A2[i][j]==(АТ *А+α*Е)

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

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

      _A2[i][j] = _One1[i][j]+_A1[i][j];

находим обратную матрицу _A2[i][j]== +=_Ar[i][j] ==(АТ *А+α*Е)-1

    this->_REV_MATR();

находим результирующий столбец _Z[i]== (АТ *А+α*Е)-1Т

    for(int i=0;i<size;i++){

     _S = 0;

     for(int j=0;j<size;j++){

        _S+=_Ar[i][j]*_B[j];

        _Z[i] = _S;

      }

     }

 

находим матрицу произведения главной матрицы и результирующего столбца _U1[i], т.е. _U1[i] == (АТ *А+α*Е)*Хα    

    for(int i=0;i<size;i++){

     _S = 0;

     for(int j=0;j<size;j++){

      _S+=_A[i][j]*_Z[j];

          _U1[i] = _S;

         }

     }

    q++;

     }

пока не выполнится условие |АТ*В- (АТ *А+α*Е)*Хα|<ε будем  продолжать цикл, увеличивая на каждом шаге k на единицу

    while(!_ALL_RIGHT());

 

}

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

При запуске программы появится окно.

9

8

7

6

5

4

3

2

1

  1.  Исходная матрица(система линейных уравнений).
  2.  Решение системы линейных уравнений.
  3.  Векторный столбец правой части уравнения Ax=b.
  4.  Порядок матрицы.
  5.  Интерактивно задаваемый диапазон элементов.
  6.  Интерактивно задаваемый диапазон определителя.
  7.  Загрузка  матрицы из файла.
  8.  Данные о программе.
  9.  Интерактивно задаваемая точность.

Вывод.

Был разработан программный продукт, реализующий метод регуляции А.Н. Тихоного для плохо обусловленных систем, учитывающий все особенности данного метода.


 

А также другие работы, которые могут Вас заинтересовать

54697. Осінь золота 237 KB
  Мета: розповсюдження поглиблення і систематизація знань учнів про пору року осінь формувати вміння працювати в групі; виховувати інтерес до літератури рідної мови та любов до природи. Обладнання: Колоски пшениці жита композиції з природніх матеріалів малюнки на тему Осінь. Осінь ведуча одягнена у однокольорове плаття на якому нашите листя клена. Пісня Школа На сцену виходять 5 учнів: 1й Заходить Осінь вклоняється тим хто в залі Непомітно зявилася осінь Все коротшає день щодобиГлянь берізки уже злотокосіІ в дубів...
54698. Антимонопольная политика государства 22.83 KB
  Рынок действует по определенным принципам, которые монополия подрывает. Поэтому борьба с монополией является одновременно защитой основных принципов рыночной экономики.
54699. Королева Осінь 37 KB
  Діти заходять під марш. Вихователь: Діти сьогодні ми з вами йдемо до країни Осінь. Під тиху музику діти йдуть на півпальцях. Діти розслаблюються.
54700. Основи безпеки життєдіяльності на уроках фізичної культури 326.5 KB
  Техніка безпеки на уроках фізичної культури 1.1 Основи техніки безпеки на уроках фізичної культури 1.2 Профілактика травматизму як основний напрямок техніки безпеки на уроках фізичної культури 1.
54701. Социально-экономическое понятие предприятия 56.5 KB
  Предприятие – это самостоятельный хозяйствующий субъект, созданный предпринимателем или объединением предпринимателей для производства продукции, выполнения работ и оказания услуг с целью удовлетворения общественных потребностей и получения прибыли
54702. Небезпека від вогню 87 KB
  Складання плану уроку: Історія вогню Що таке пожежа Причини виникнення пожежі Невідкладні дії під час пожежі Як поводитися під час пожежі ІІ. Вчитель: Ця група нагадала нам телефон пожежної служби яка прибуває на допомогу під час пожежі але пожежники кажуть: Пожежі краще запобігти ніж її гасити. Тому передаємо слово групі Юні пожежні які нагадають вам причини виникнення пожежі. Вона показує нам основні джерела займання та причини виникнення пожежі.
54703. Основы генетики, урок 71.5 KB
  Урок: освоение нового материала. Вильгельм Швебель немецкий ученый и публицист Слайд 2 Ход урока: Организационный момент Актуализация знаний. Слайд 3 Генетика рассматривает носители наследственности хромосомы и гены. Учащиеся записывают в тетради Слайд 4 Наследственностью обладают все живые организмы.
54705. Сценарій Останнього уроку 123 KB
  Шановні гості Ми сьогодні зібралися в цій залі щоб вирядити на широку дорогу дорослого життя наших випускників.її життя починається світанком душі дитинством. Це справді найкраща пора життя людини і проходить вона у школі. Дороги дороги Тернисті і рівні Із пристані школа ведуть у життя.