99703

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

Курсовая

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

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

Русский

2016-10-08

117.26 KB

4 чел.

Содержание

  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.  Интерактивно задаваемая точность.

Вывод.

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


 

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

41099. Управление заемным капиталом 1.25 MB
  Обеспечение своевременных расчетов по полученным кредитам На второй стадии анализа определяются основные формы привлечения заемных средств анализируются в динамике удельный вес сформированных финансового кредита товарного кредита и текущих обязательств по расчетам в общей сумме заемных средств используемых предприятием. Эти формы дифференцируются в разрезе финансового кредита; товарного коммерческого кредита; прочих форм. К числу важнейших из этих условий относятся; а срок предоставления кредита; б ставка процента за кредит;...
41101. ОСНОВЫ МЕЖДУНАРОДНЫХ ВАЛЮТНЫХ, РАСЧЕТНЫХ И КРЕДИТНЫХ ОТНОШЕНИЙ 370.33 KB
  С середины 70х годов используется паритет на базе валютной корзины; режим валютного курса фиксированный плавающий в определенных пределах; наличие или отсутствие валютных ограничений; регулирование международной валютной ликвидности страны; официальные золотые и валютные резервы стран; счета СДР ЭКЮ резервная позиция в МВФ; регламентация использования международных кредитных средств обращения и форм международных расчетов; режимы валютного рынка и рынка золота;статус национальных органов регулирующих валютные...
41102. Центральные и коммерческие банки, их функции и формы организации 292.94 KB
  Введение Термин банк происходит от итальянского слова банко что означает лавка скамья или конторка за которой менялы оказывали свои услуги. Во многих источниках дошедших до нас можно встретить данные о вавилонских банкирах принимавших процентные вклады и выдававших ссуды под письменные обязательства и под залог различных ценностей. Вавилонский банк принимал вклады платил по ним проценты выдавал ссуды и даже выпускал банковские билеты.
41103. Компьютерные мониторы на основе электронно-лучевой трубки 839 KB
  Сквозь металлическую маску или решетку они попадают на внутреннюю поверхность стеклянного экрана монитора которая покрыта разноцветными люминофорными точками. Поток электронов луч может отклоняться в вертикальной и горизонтальной плоскости что обеспечивает последовательное попадание его на все поле экрана. Чтобы электроны беспрепятственно достигали экрана из трубки откачивается воздух а между пушками и экраном создаётся высокое электрическое напряжение ускоряющее электроны.Это сделано для того чтобы электронный луч в центре экрана и...
41104. Цифровая бумага 1.87 MB
  Многоцветная полихромная электронная бумага Электронная бумага EDP В отличие от традиционных жидкокристаллических плоских дисплеев в которых используется просвет матрицы для формирования изображения электронная бумага формирует изображение в отраженном свете как обычная бумага и может показывать текст и графику неопределенно долго не потребляя при этом электричество и позволяя изменять изображение в дальнейшем.
41105. ПРОЦЕССЫ САМОТЕСТИРОВАНИЯ КОМПЬЮТЕРНОЙ СИСТЕМЫ ПРИ ВКЛЮЧЕНИИ 657.5 KB
  Блок регистров общего назначения определяет вычислительные ресурсы микропроцессора и содержит регистры для временного хранения данных и команд. Данные и команды передаются по шине данных а шина управления состоит из линий для передачи управляющих сигналов состояния памяти и периферийных устройств. С точки зрения структуры микропроцессора для пользователя присутствуют следующие параметры: архитектура адресное пространство памяти достижимое процессором разрядность шины данных быстродействие. Передача данных в режиме прямого доступа к...
41107. Виды обеспечения АСУ. Техническое обеспечение АСУ 32.5 KB
  Организационное обеспечение АСУ представляет собой совокупность средств и методов, предназначенных для проведения технико-экономического анализа существующей системы управления, выбора и постановки задач автоматизации организационного управления предприятием, организации производства и управления в условиях АСУ.