99703

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

Курсовая

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

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

Русский

2016-10-08

117.26 KB

2 чел.

Содержание

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

Вывод.

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


 

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

35907. Классы неорганических соединений 76 KB
  Кислотные и основные оксиды реагируют между собой образуя соли: CO2K2O = K2CO3; SiO2CO = CSiO3. Кислоты реагируют с основными оксидами с образованием соли и воды: 2HClMgO = MgCl2H2O. Кислоты взаимодействуют с основаниями с образованием соли и воды реакция нейтрализации: H2SO4CuOH2 = CuSO42H2O. Кислоты реагируют с амфотерными гидроксидами оксидами с образованием соли и воды: 2HClZnOH2 = ZnCl22H2O.
35908. Технология Токайских вин. Изменения происходящие в ягоде при поражении ее плесневым грибом «Botrytis cinerea» (серая гниль). Их влияние на ход дальнейших технологических операций 74.5 KB
  Токай венгерское вино готовится из сортов винограда Фурминт главным образом Гарс Левелю и Мускат. Особенностью технологии токайских вин Венгрии является использование наряду с виноградом зрелым или слегка перезрелым также увяленных и заизюмленных ягод пораженных грибом Ботритис Цинереа. Считают что благодаря воздействию этого гриба в токайских винах образуются специфические букет и вкус.
35909. Местное самоуправление в странах европейского союза: особенности правового регулирования и специфика функционирования 74.5 KB
  местное самоуправление в странах европейского союза: особенности правового регулирования и специфика функционирования Регулируется Европейской хартией местного самоуправления которая принята Советом Европы 15 октября 1985 года . Прежде всего полномочия не могут предоставляться органам местного самоуправления на временной основе они должны основываться на национальном законодательстве. Таким образом компетенция органов местного самоуправления должна определяться по мере возможности Конституцией и законодательством. Органы местного...
35910. Поток вектора магнитной индукции. Теорема Гауса для вектора магнитной индукции 74.05 KB
  Контур с током в магнитном поле. Взаимодействие токов осуществляется через поле которое называется магнитным. Следовательно движущиеся заряды токи изменяют свойства окружающего их пространства создают в нем магнитное поле. Это поле проявляется в том что на движущиеся в нем заряды токи действуют силы.
35911. Информационные системы. ИТ 74 KB
  Термин информационные системы Информационная система ИС – это организационноупорядоченная взаимосвязанная совокупность средств и методов ИТ а также используемых для хранения обработки и выдачи информации в интересах достижения поставленной цели. Такое понимание информационной системы предполагает использование в качестве основного технического средства переработки информации ЭВМ и средств связи реализующих информационные процессы и выдачу информации необходимой в процессе принятия решений задач из любой области. Таким образом ИТ...
35913. Экология как наука 73.5 KB
  Все живые и неживые объекты окружающие растения животные и другие организмы непосредственно взаимодействующие с ними называются средой обитания. Световой режим оказывает прямое влияние в первую очередь на растения. По отношению к освещенности выделяют следующие экологические группы растений: – гелиофиты – светолюбивые растения растения открытых пространств постоянно хорошо освещаемых местообитаний. – сциофиты – тенелюбивые растения которые плохо переносят интенсивное освещение растения нижних ярусов тенистых лесов.
35914. Государственное регулирование экономики: сущность, средства, объекты. Институционализм 75.5 KB
  Какие же проблемы не может решить рынок 1. А поскольку сегодня практически во всех странах мира преобладает рынок несовершенной конкуренции то регулирующая роль государства значительно возрастает. Оно должно обеспечивать условия для эффективного развития производства; перераспределять доходы хозяйствующих субъектов и населения в целом; оказывать поддержку малообеспеченным; регулировать рынок рабочей силы; смягчать цикличность развития экономики; проводить антиинфляционную и антикризисную политику; способствовать развитию конкуренции...
35915. Советский плакат 71.12 KB
  Политический плакат Основная статья: Советский плакат Виды искусства способные жить на улицах в первые годы после революции играли важнейшую роль в формировании общественного и эстетического сознания революционного народа. Он оказался самым мобильным и оперативным видом искусства[4]. В его плакатах заметна символикоаллегорическая тенденция советского искусства той поры.