Реферат: Методы решения некорректно поставленных задач

W[ z ]<=d, компактно на F.

3.2.2. Итак, рассмотрим произвольную систему линейных алгебраических уравнений (короче — СЛАУ)

Аz =u,                                            (3; 2,2)

в которой z и u—векторы, z=(z1, z2, ...,zn) ÎRn, и=(u1,u2, ...,un) ÎRm, А—матрица с элементами aij, А = {aij}, где j =1, 2, ..., n; i= 1, 2, ..., т, и число п не обязано быть равным числу т.

Эта система может быть однозначно разрешимой, вы­рожденной (и иметь бесконечно много решений) и не­разрешимой.

Псевдорешением системы (3; 2,2) называют вектор z, минимизирующий невязку || Az – u || на всем пространстве Rn. Система (3; 2,2) может иметь не одно псевдоре­шение. Пусть FA есть совокупность всех ее псевдореше­ний и z1 — некоторый фиксированный вектор из Rn, оп­ределяемый обычно постановкой задачи.

Нормальным относительно вектора z1 решением си­стемы (3;2,2) будем называть псевдорешение z0 с ми­нимальной нормой || z – z1 ||, т. е. такое, что

|| z0 – z1 || =

Здесь  . В дальнейшем для простоты записи будем считать z1= 0 и нормальное относительно вектора z1=0 решение называть просто нормальным ре­шением.

Для любой системы вида (3; 2,2)  нормальное решение существует и единст­венно.

Замечание 1. Нормальное решение z° системы (3;2,2) можно определить также как псевдорешение, минимизирующее заданную положительно определенную квадратичную форму относительно координат вектора z—z1. Все излагаемые ниже результаты остаются при этом справедливыми.

Замечание 2. Пусть ранг матрицы А вырожден­ной системы (3; 2,1) равен r < n и zr+1,zr+2, … , zn— базис линейного пространства NA, состоящего из элемен­тов z, для которых Аz=0, NA = {z; Аz= 0}. Решение z° системы (3; 2,1), удовлетворяющее n—r условиям ортогональности

                                 (z0 – z1, zS)= 0,   S= r + 1, r + 2, .. ,n,                      (3; 2,3)

определяется однозначно и совпадает с нормальным ре­шением.

3.2.3. Нетрудно видеть, что задача нахождения нормаль­ного решения системы (3; 2,2) является некорректно поставленной. В самом деле, пусть А — симметричная матрица. Если она невырожденная, то ортогональным преобразованием

z = Vz*, u = Vu*

ее можно привести к диагональному виду и преобразо­ванная система будет иметь вид

lizi*=ui* , i= 1, 2,. .., п,

где li — собственные значения матрицы А.

Если симметричная матрица А — невырожденная и имеет ранг r, то n – r  ее собственных значений равны нулю. Пусть

li¹0 для i=1, 2, ..., r;

и

li=0 для i=r+1,r+2, …, n.

Полагаем, что система (3; 2,2) разрешима. При этом ui*= 0 для i =r + 1, ..., n.

Пусть «исходные данные» системы (А и и) заданы с погрешностью, т. е. вместо А и и заданы их прибли­жения А и u:

 || A – A ||<=h,  ||u – u||<=d . При этом           

                                                        (3;2,4)

Пусть li — собственные значения матрицы А. Извест­но, что они непрерывно зависят от А в норме (3; 2,4). Следовательно, собственные значения lr+1 , lr+2, …,ln могут быть сколь угодно малыми при достаточно малом h.

Если они не равны нулю, то

             zi*=.

Таким образом, найдутся возмущения системы в пре­делах любой достаточно малой погрешности А и и, для которых некоторые zi* будут принимать любые наперед заданные значения. Это означает, что задача нахожде­ния нормального решения системы (3; 2,2) является неустойчивой.

Ниже дается описание  метода нахождения нормального   решения системы (3; 2,2), устойчивого к малым (в норме (3; 2,4)) возму­щениям правой части и , основанного на методе регуляризации.

3.3. Метод регуляризации нахождения нормального решения

3.3.1. Пусть z° есть нормальное решение системы

Аz = и.                   (3; 3,1)

Для простоты будем полагать, что прибли­женной может быть лишь правая часть, а оператор (матрица) А — точный.

Итак, пусть вместо и мы имеем вектор и, || и — и || <=d ; т. е. вместо системы (3;3,1) имеем систему

(3; 3,2)

Аz = u.

Требуется найти приближение zd к нормальному реше­нию системы (3;3,1), т. е. к вектору z° такое, что zd àz° при d à0. Отметим, что векторы u и u (один из них или оба) могут не удовлетворять классическому ус­ловию разрешимости.

Поскольку система (3; 3,1) может быть неразреши­мой, то inf ||Az-u|| = m >=0, где inf берется по всем векторам z Î Rn.

Естественно искать приближения zd  в классе Qd  век­торов z, сопоставимых по точности с исходными данны­ми, т. е. таких, что || Az – u ||<=m+d. Но поскольку вместо вектора u мы имеем вектор u, то мы можем найти лишь

   m=inf || Az – u ||.

 zÎ Rn

Отметим, что из очевидных неравенств

||Az – u ||<=||Az – u || + || u – u || ,

                                                            ||Az – u ||<= || Az – u || + ||u – u ||

следуют оценки m<=m+d, m<=m+d, приводящие к не­равенству | m — m | <=d. Поэтому будем искать прибли­жение zd к нормальному решению z° в классе Qd векто­ров z, для которых || Аz — и || <=m +2d. Отметим, что если имеется информация о разрешимости системы (3;3,1), то m = 0 и в качестве класса Qd можно брать класс векторов z, для которых || Аz — и|| <= d. Класс Qd есть класс формально возможных приближенных решений.

Но нельзя в качестве zd брать произвольный вектор из класса Qd, так как такое «приближение» будет неустойчивым к малым изменениям правой части уравнения (3;3,2). Необходим принцип отбора. Он естественным образом вытекает из постановки задачи. В самом деле согласно определению нормального решения искомое ре­шение z° должно быть псевдорешением с минимальной нормой. Поэтому в качестве приближения к z° естествен­но брать вектор zd  из Qd, минимизирующий функционал

W[ z ] = ||z||2 на множестве Qd .

Таким образом, задача сводится к минимизации функ­ционала W[ z ] = ||z||2  на множестве Qd векторов z, для которых выполняется условие || Аz — u || <=m +2d.

3.3.2. Пусть zd  — вектор из Qd, на котором функционал ||z||2 достигает минимума на множестве Qd. Его можно рассматривать как результат применения к правой части u уравнения (3; 3,2) некоторого оператора R1(u, d), зависящего от параметра d. Справедлива

Теорема 1. Оператор R1(u, d) обладает следующи­ми свойствами:

1) он определен для всякого uÎRm и любого d > 0;

2) при d à0  zd== R1(u, d)  стремится к нормальному решению z° уравнения Аz=u, т. е. он является регуляризирующим для уравнения Аz=u .

3.3.3. Пусть  zd  — вектор, на котором функционал W[ z ] = ||z||2   достигает минимума на множестве Qd. Легко ви­деть из наглядных геометрических представлений, что вектор zd  лежит на границе множества Qd, т.е. ||Azd - u ||=m +2d =d1.

 Это следует непосредственно также из того, что функционал W[ z ] = ||z||2  является сстабилизирующим и квазимонотонным. Стабилизирующий функционал W[ z ] называется квазимонотонным , если каков бы ни был элемент z из F1 , не принадлежащий множеству M0 , в любой его окрестности найдется элемент z1 из F1, для которого W[ z1 ]< W[ z ], т.е. если функционал не имеет локальных минимумов на множестве F1\ M0.

Задачу нахождения вектора zd   можно поставить так: среди векторов z, удовлетворяющих условию ||Az – u ||=m +2d  , найти вектор zd  с минимальной нормой, т. е. минимизирующий функционал W[ z ]=||z||2.

Последнюю задачу можно решать методом Лагранжа, т. е. в качестве zd брать вектор za, минимизирующий функционал

Мa [z, u] = ||Az - u ||2+ a||z||2,         a>0,

с параметром a, определяемым по невязке, т. е. из ус­ловия ||Аza— u||=d1.  При этом параметр  a определяется однозначно .

3.3.4. Поскольку Мa [z, u]  — квадратичный функционал, то для любых u ÎRm и  a> 0 существует лишь один минимизирующий его вектор za. В самом деле, допустим,

что существуют два вектора za и za, минимизирующие его. Рассмотрим векторы z, расположенные на прямой (пространства Rn), соединяющей za и za:

z = za + b( za - za).

Функционал Мa [z, u]  на элементах этой прямой есть неотрицательная квадратичная функция от b. Следова­тельно, она не может достигать наименьшего значения при двух различных значениях b: b = 0 (z = za) и b=1 (z = za).

Компоненты zja вектора za  являются решением си­стемы линейных алгебраических уравнений

получающихся из условий минимума функционала Мa [z, u]:

Здесь

Компоненты zja  могут быть определены и с помощью какого-нибудь другого алгоритма минимизации            функцио­нала Мa [z, u].

Вектор za   можно рассматривать как результат приме­нения к u некоторого оператора za=R(u, a), завися­щего от параметра a.

Покажем, что оператор R0(u, a) является регуляризирующим для системы (3;3,1), т. е. обладает свойствами 1) и 2) определения 2 (см. 3.1.2.). В п. 3.3.2. было сказано, что он определен для всяких u ÎRm и  a > О и, следовательно, обладает свойством 1). Теперь пока­жем справедливость свойства 2), т. е. существование таких функций a=a(d) , что векторы za(d) = R0(u, a(d)) сходятся к нормальному решению z° системы (3; 3,1) при dà0. Это непосредственно следует из приводимой ниже теоремы 2.

Теорема 2( Тихонова). Пусть z° есть нормальное решение си­стемы Az= u и вместо вектора u мы имеем вектор u  такой, что ||u—u||<=d. Пусть, далее, b1(d) и b2(d) — какие-либо непрерывные на [0, d2] и положительные на (0, d2] функции, монотонно стремящиеся к нулю при dà 0 и такие, что

                                     

Тогда для любой .положительной на (0, d2] функции a=a(d) , удовлетворяющей условиям

векторы za(d) = R0(u, a(d))  сходятся к нормальному ре­шению z0 системы Az = u при dà0, т. е.

Примечание. Доказательства теорем в данном разделе опущены, т.к. основной теоретической частью работы является раздел «Метод Подбора. Квазирешения». Метод Тихонова описан из-за использования его в численном эксперименте.

ЗАКЛЮЧЕНИЕ

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

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

z=A-1u

Получим    z1=316

                    z2=-990

                    z3=832

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

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

                   z1=348

                   z2=-1090

                   z3=916.

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

            Будем искать решение методом Тихонова. В теоретической части было показано, что целесообразно использовать регуляризирующий оператор следующего вида:  (aE + ATA)za=ATud , где E – единичная матрица, za -- приближенное нормальное решение, AT – транспонированная исходная матрица, a -- параметр регуляризации,

ud --  правая часть, заданная неточно. Эту задачу можно решать стандартными методами, задав предварительно функцию a=a(d) , удовлетворяющую условиям теоремы Тихонова. В моем примере это функция  a(d)=d/4d. Далее будем решать регуляризованную задачу с точностью e=0.001 ,последовательно изменяя значения a.

В качестве контр-примера можно подставить в программу любую функцию a(d) , не удовлетворяющую условиям теоремы Тихонова. Любая положительная функция монотонно возрастающая, не обладающая свойством a(d)à0 при dà0, не будет минимизировать невязку.

Текст программы приведен в приложении 1. Полная распечатка результатов приведена в приложении 2. Здесь же представлены окончательные значения на выходе из программы.

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

Z(1)= 3.47834819174013E+0002

Z(2)=-1.08948394975175E+0003

Z(3)= 9.15566443137791E+0002

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

U1(1)= 9.99997717012495E-0001

U1(2)= 1.00000741970775E+0000

U1(3)= 1.09948402394883E+0000

Значение параметра регуляризации:

 2.61934474110603E-0010

ПРИЛОЖЕНИЯ

Приложение 1.

     Текст программы для реализации метода Тихонова на языке PASCAL

Uses CRT;

type

 real=extended;

const

  matrixA: array[1..3,1..3] of real  = ((-19/20,1/5,  3/5),

                                      (-1    ,0.1,  0.5),

                                      (-0.01 ,0  ,1/200));

  One: array [1..3,1..3] of real = ((1,0,0),

                                    (0,1,0),

                                    (0,0,1));

  U:array[1..3] of real = (1,1,1.1);

  var

   i,j,k,q:byte;

   A,At,A1,A2,Ar,One1:array[1..3,1..3] of real;

   delta,Det,S,alpha:real;

   B,Z,U1:array[1..3] of real;

   f:text;

   Procedure TransA;

    begin

     for i:=1 to 3 do

      for j:=1 to 3 do

         At[i,j]:=A[j,i]

    end;

   Function Koef(par1,par2:byte):real;

    var

     Sum:byte;

     Tmp:real;

      begin

       Sum:=par1+par2;

       Tmp:=1;

       for k:=1 to sum do

        Tmp:=Tmp*(-1);

       Koef:=Tmp;

      end;

   Function AlAdd(par1,par2:byte):real;

    type

     element=record

             value:real;

             flag:boolean;

             end;

    var

     BB:array[1..2,1..2] of real;

     AA:array[1..3,1..3] of element;

     k,v,w:byte;

     N:array[1..4] of real;

     P1:real;

    begin

     for v:=1 to 3 do

      for w:=1 to 3 do begin

       AA[v,w].value:=A2[v,w];

       AA[v,w].flag:=true

      end;

     for v:=1 to 3 do AA[par1,v].flag:=false;

     for v:=1 to 3 do AA[v,par2].flag:=false;

     { for v:=1 to 3 do begin

      for w:=1 to 3 do write(AA[i,j].value:2:3,' ');

      writeln

     end; }

     k:=1;

     for v:=1 to 3 do

      for w:=1 to 3 do

       begin

       if AA[v,w].flag then

        begin

         N[k]:=AA[v,w].value;

      {   writeln(N[k]);}

         k:=k+1

        end;

       end;

     BB[1,1]:=N[1];  BB[1,2]:=N[2];

     BB[2,1]:=N[3];  BB[2,2]:=N[4];

{     writeln('alg dop',par1,par2,' ',BB[1,1]*BB[2,2]-BB[1,2]*BB[2,1]);}

     AlAdd:=BB[1,1]*BB[2,2]-BB[1,2]*BB[2,1];

    end;

    Function DetCount:real;

     var

       S1:real;

       z:byte;

     begin

      S1:=0;

      for z:=1 to 3 do S1:=S1+A2[1,z]*Koef(1,z)*AlAdd(1,z);

      DetCount:=S1;

     end;

    Procedure RevMatr;

     begin

       for i:=1 to 3 do

        for j:=1 to 3 do

         Ar[j,i]:=Koef(i,j)*AlAdd(i,j)/DetCount;

{       for i:=1 to 3 do begin

        for j:=1 to 3 do write(Ar[i,j],' ');

        writeln;

       end;}

     end;

     Function AllRight:boolean;

      begin

       writeln(f,'­Ґўп§Є  Ї® 1-¬г н«-вг',(abs(U[1]-U1[1])));

       writeln(f,'­Ґўп§Є  Ї® 2-¬г н«-вг',(abs(U[2]-U1[2])));

       writeln(f,'­Ґўп§Є  Ї® 3-¬г н«-вг',(abs(U[3]-U1[3])));

       writeln(F);

       if (abs(U[1]-U1[1])<0.001) and (abs(U[2]-U1[2])<0.001) and

          (abs(U[3]-U1[3])<0.001) then AllRight:=true

          else AllRight:=false

      end;

     Function Pow(par1:real;par2:byte):real;

      var

       S2:real;

       z:byte;

      begin

       S2:=1;

       if par2=0 then begin

        Pow:=1;

        exit

       end

         else

           for z:=1 to par2 do S2:=S2*par1;

       Pow:=S2;

      end;

      BEGIN

       clrscr;

       Assign(f,'c:\tikh.txt');

       Rewrite(f);

       for i:=1 to 3 do

        for j:=1 to 3 do

         A[i,j]:=matrixA[i,j];

       TransA;

       Det:=0.000125;

       {----------------------------}

       for i:=1 to 3 do begin

        S:=0;

        for j:=1 to 3 do begin

          S:=S+At[i,j]*U[j];

          B[i]:=S

         end;

       end;

       {----------------------------}

       for i:=1 to 3 do

        for j:=1 to 3 do

         begin

          S:=0;

           for k:=1 to 3 do begin

            S:=S+At[i,k]*A[k,j];

            A1[i,j]:=S

           end

         end;

       {-----------------------------}

       q:=1;

      repeat

       alpha:=q/pow(4,q);

       for i:=1 to 3 do

        for j:=1 to 3 do

         One1[i,j]:=One[i,j]*alpha;

       for i:=1 to 3 do

        for j:=1 to 3 do

         A2[i,j]:=One1[i,j]+A1[i,j];

       RevMatr;

       {------------------------------}

       for i:=1 to 3 do begin

        S:=0;

        for j:=1 to 3 do begin

          S:=S+Ar[i,j]*B[j];

          Z[i]:=S

         end;

       end;

       for i:=1 to 3 do begin

        S:=0;

         for j:=1 to 3 do begin

          S:=S+A[i,j]*Z[j];

          U1[i]:=S

         end

       end;

       q:=q+1;

       until AllRight;

       {------------------------------}

       clrscr;

       writeln('ЏаЁЎ«Ё¦Ґ­ЁҐ Є ­®а¬ «м­®¬г аҐиҐ­Ёо');

       for i:=1 to 3 do writeln('Z(',i,')=',z[i]);

       writeln;

       writeln('‡­ зҐ­ЁҐ Їа ў®© з бвЁ ЇаЁ Ї®¤бв ­®ўЄҐ ЇаЁЎ«. аҐиҐ­Ёп');

       for i:=1 to 3 do writeln('U1(',i,')=',U1[i]);

       writeln;

       writeln('‡­ зҐ­ЁҐ Ї а ¬Ґва  аҐЈг«паЁ§ жЁЁ:');

       writeln(alpha);

       Close(f);

       readln;

      END.

Приложение 2.

   Распечатка результатов пересчета на каждом шаге

невязка по 1-му эл-ту 7.75620788018006E-0002

невязка по 2-му эл-ту 9.12970302562861E-0002

невязка по 3-му эл-ту 1.09101153877771E+0000

невязка по 1-му эл-ту 3.51667654246499E-0002

невязка по 2-му эл-ту 4.81631787337596E-0002

невязка по 3-му эл-ту 1.09057642915500E+0000

невязка по 1-му эл-ту 8.14255746519741E-0003

невязка по 2-му эл-ту 1.75271999674588E-0002

невязка по 3-му эл-ту 1.09024740493812E+0000

невязка по 1-му эл-ту 1.64128226088452E-0004

невязка по 2-му эл-ту 1.40420815653456E-0003

невязка по 3-му эл-ту 1.09002512985506E+0000

невязка по 1-му эл-ту 1.09651876415789E-0003

невязка по 2-му эл-ту 8.01044623892439E-0003

невязка по 3-му эл-ту 1.08980075500722E+0000

невязка по 1-му эл-ту 3.24092274239579E-0003

невязка по 2-му эл-ту 1.28969442769472E-0002

невязка по 3-му эл-ту 1.08943309314635E+0000

невязка по 1-му эл-ту 4.29878415191160E-0003

невязка по 2-му эл-ту 1.47864580098997E-0002

невязка по 3-му эл-ту 1.08840358157784E+0000

невязка по 1-му эл-ту 4.64764022304719E-0003

невязка по 2-му эл-ту 1.53489294761093E-0002

невязка по 3-му эл-ту 1.08488736141985E+0000

невязка по 1-му эл-ту 4.70263264899617E-0003

невязка по 2-му эл-ту 1.53524096326819E-0002

невязка по 3-му эл-ту 1.07252416252061E+0000

невязка по 1-му эл-ту 4.54618391386039E-0003

невязка по 2-му эл-ту 1.47935415193105E-0002

невязка по 3-му эл-ту 1.03007092553528E+0000

невязка по 1-му эл-ту 3.97950585276394E-0003

невязка по 2-му эл-ту 1.29378307693635E-0002

невязка по 3-му эл-ту 9.00028069734717E-0001

невязка по 1-му эл-ту 2.71895340473448E-0003

невязка по 2-му эл-ту 8.83742514077426E-0003

невязка по 3-му эл-ту 6.14624514462952E-0001

невязка по 1-му эл-ту 1.25089904346179E-0003

невязка по 2-му эл-ту 4.06552487723671E-0003

невязка по 3-му эл-ту 2.82729125073012E-0001

невязка по 1-му эл-ту 4.15581257891512E-0004

невязка по 2-му эл-ту 1.35064829843828E-0003

невязка по 3-му эл-ту 9.39264706989556E-0002

невязка по 1-му эл-ту 1.18814900667952E-0004

невязка по 2-му эл-ту 3.86149131520602E-0004

невязка по 3-му эл-ту 2.68533566153482E-0002

невязка по 1-му эл-ту 3.22671215741144E-0005

невязка по 2-му эл-ту 1.04868192738639E-0004

невязка по 3-му эл-ту 7.29267248287954E-0003

невязка по 1-му эл-ту 8.61328853146714E-0006

невязка по 2-му эл-ту 2.79931897352870E-0005

невязка по 3-му эл-ту 1.94668264668650E-0003

невязка по 1-му эл-ту 2.28298750498679E-0006

невязка по 2-му эл-ту 7.41970775380851E-0006

невязка по 3-му эл-ту 5.15976051172231E-0004

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

Z(1)= 3.47834819174013E+0002

Z(2)=-1.08948394975175E+0003

Z(3)= 9.15566443137791E+0002

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

U1(1)= 9.99997717012495E-0001

U1(2)= 1.00000741970775E+0000

U1(3)= 1.09948402394883E+0000

Значение параметра регуляризации:

 2.61934474110603E-0010

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1.А.Н.ТИХОНОВ, В.Я.АРСЕНИН  «МЕТОДЫ РЕШЕНИЯ НЕКОРРЕКТНЫХ ЗАДАЧ» – МОСКВА «НАУКА» 1979.

2.Г.И.МАРЧУК  «МЕТОДЫ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ» – МОСКВА «НАУКА» 1977.

3.Л.И.ГОЛОВИНА «ЛИНЕЙНАЯ АЛГЕБРА И НЕКОТОРЫЕ ЕЕ ПРИЛОЖЕНИЯ» – МОСКВА «НАУКА»  1975.

4.В.И.РАКИТИН, В.Е.ПЕРВУШИН  «ПРАКТИЧЕСКОЕ РУКОВОДСТВО ПО МЕТОДАМ ВЫЧИСЛЕНИЙ» – МОСКВА «ВЫСШАЯ ШКОЛА» 1998.

5.В.В.ФАРОНОВ «ПРОГРАММИРОВАНИЕ НА ПЕРСОНАЛЬНЫХ ЭВМ В СРЕДЕ TURBO PASCAL»  -- ИЗДАТЕЛЬСТВО МГТУ  1990.