Форум на Исходниках.RU https://forum.sources.ru Форум на Исходниках.RU Форум на Исходниках.RU https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=759610 Sun, 26 Jun 2005 12:35:08 +0000 Численные методы https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=759610 Вычисление приближённого значения функций с помощью ряда Тейлора (Разложение функций в степенные ряды) (применяется также для нахождения пределов вообще)

К примеру, sin = x -(x3/3!)+(x5/5!)- ...
Здесь описано как вычислять Sin, Cos, ArcTg, Exp. Другие функции пишутся аналогично.
Значение функции вычисляется в процедуре FunctionValue. Только подставьте нужное название функции.


Для того, чтобы вычислить приближённое значение по формуле, мы должны взять бесконечно большое число, чтобы приблизиться к истинному значению функции. Т.е. при n->infinity (бесконечность), получаем истинное значение функции в данной точке.
На практике мы этого не станем делать, так как в Паскале имеются типы с ограниченным размером, такие как: Real, Single, Double, Extended. Что же теперь делать?

Введём понятие точность вычисления и назовём эту переменную eps (epsilon). Точность вычисления может быть разная: 0.05, 0.001 и намного большая, 10-8. В большинстве случаев хватает точности 3 знака после нуля, поэтому возьмём eps=0.001, или в краткой форме записи для Паскаля, 1e-3, что означает 1*10-3. Типа Single для этой точности вполне хватит, на нём и остановимся.
Если истинное значение функции в заданной точке равно, допустим, 1.2328, а вычисленное значение равно 1.2323, то при нашей точности вычисления разница по модулю между ними составляет 0.0005<eps, то есть ответ нас устраивает.

Для вычисления бесконечного ряда нужно применить цикл, в котором мы каждый раз будем проверять, насколько точно вычислено требуемое значение X, то есть проверим:|Xnew - Xold| < eps ?
где Xnew - это новое вычисленное значение,
Xold - предыдущее вычисленное значение
Как только в ходе вычисления внутри цикла мы выясним, что достигнута требуемая точность (eps), выходим из него. Дальше вычислять уже нет необходимости.

Далее. Как вычислять эти значения? Смотрим на общую формулу:
sin (x) = Sum [(-1)n * x2n+1 / (2n+1)!] , где n->infinity.
Из чего состоит данная формула?
Из вычисления степени числа (функция Power), факториала числа (функция Factorial) и вычисления знака (функция Sign).
Для повторения вычисления этих функций смотрите прикреплённые вверху ссылки.

Теперь надо собрать все эти функции в одном выражении, которое будет вычисляться в отдельной функции.
     { вычисление значения элемента с индексом [I]i[/I] ряда Тейлора для синуса}
    function Sine (x: single; i: integer): single;
    var k: integer;
    begin
      K := 2 * i + 1;
      Sine := ( Sign (i) * Power (x, k) / factorial (k) );
    end;

Что осталось? Объединить всё написанное в одной функции вычисления приближённого значения FunctionValue.
    function FunctionValue (x: single): single;
На входе задаётся точка (число), в которой нужно вычислить значение функции.

А теперь программа:

    const eps: single = 1e-3; { Epsilon - необходимая точность, 1*10^(-3) }
     
    function Factorial (N: word): single; { N! }
    var
      f: single;
    begin
      Factorial := 1.0;
      if n = 0 then exit;
     
      f := 1.0;
      for n := 1 to n do
        f := f * n;
     
      Factorial := f;
    end;
     
    function Sign (p: integer): single; { (-1)^p }
    begin
         if Not Odd(P) then
           Sign := 1
         else
           Sign := -1;
    end;
     
    function Power (base: single; N: integer): single; { степень N по основанию base }
    var k: integer;
        P: single;
    begin
         P := 1.0;
         for k := 1 to N do
           P := P * base;
     
         Power := P;
    end;
     
    {--------------------------------------------------------------------}
     
    function Sine (x: single; i: integer): single;
    var k: integer;
    begin
      K := 2 * i + 1;
      Sine := ( Sign (i) * Power (x, k) / factorial (k) );
    end;
     
    function CoSine (x: single; i: integer): single;
    var k: integer;
    begin
      K := 2 * i;
      CoSine := ( Sign (i) * Power (x, k) / factorial (k) );
    end;
     
    function ArcTg (x: single; i: integer): single;
    var k: integer;
    begin
      K := 2 * i + 1;
      ArcTg := ( Sign (i) * Power (x, k) / k );
    end;
     
    function Exp (x: single; i: integer): single;
    begin
      Exp := Power (x, i) / factorial (i);
    end;
     
    {--------------------------------------------------------------------}
     
    function FunctionValue (x: single): single;
    var sum,old: single;
        index: integer;
    begin
         if x <= 0.0 then
         begin
              FunctionValue := 0.0;
              exit;
         end;
         index := 0;
         sum := 0.0;
         repeat { повторять }
               old := sum; { old - прежнее вычисленное значение функции }
               sum := sum + ArcTg (x, index);   { Sin, Cos, Exp }
               index := index + 1;
         until abs (sum - old) < eps; { пока не достигнута необходимая точность }
         FunctionValue := sum
    end;
     
    begin
         writeln ( 'ArcTg (x) =', FunctionValue ( 1.0 / sqrt (3.0) ) : 8 : 3 );  { ArcTg }
    {     writeln ( 'Cosinus (x) =', FunctionValue ( pi / 6.0 ) : 8 : 3 );  { CoSine }
         readln;
    end.
]]> Romtek Pascal: Математика https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=751240 Thu, 16 Jun 2005 21:36:31 +0000 Численные методы https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=751240 Методы решения дифференциальных уравнений.

Данная статья посвящена численным методам решения дифференциальных уравнений.
Мы будем рассматривать методы решения одного обыкновенного дифференциального
уравнения первого порядка с одним начальным условием:

        y' = f(x,y)                         (1)
        y(x[0]) = y[0]                          (2)


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

В основном существуют два широких класса методов:
  • Одноступенчатые методы, в которых используется только информация о самой
    кривой в одной точке и не производятся итерации. Один из методов - решение с
    помощью рядов Тейлора, но на практике он не слишком удобен для использования.
    Практически удобные методы этого класса - методы Рунге-Кутта. Эти методы
    прямые (без итераций), но требуют многократных повторных вычислений функции.
    При использовании данных методов трудно оценивать допускаемую ошибку.
  • Многоступенчатые методы, в которых следующую точку кривой можно найти, не
    производя так много повторных вычислений, но для достижения достаточной
    точности требуются итерации. Большинство методов этого класса называются
    методами прогноза и коррекции (метод Адамса-Бошфора). некоторые трудности,
    связанные с использованием итерационной процедуры и с использованием
    нескольких начальных точек уравновешиваются тем фактом, что оценку ошибки
    при использовании этого метода легко получить в качестве побочного продукта
    вычислений.

Как и во многих других случаях, эти два класса методов придется сочетать
разумным образом, учитывая их достоинства и недостатки.

1. Методы Рунге-Кутта.

Методы Рунге-Кутта обладают следующими отличительными свойствами:
  • Эти методы одноступенчатые: чтобы найти y[m+1], нужна информация только о
    предыдущей точке x[m],y[m].
  • Они согласуются с рядом Тейлора вплоть до членов порядка h^p,
    где p - различна для разных методов и называется порядком метода.
  • Они не требуют вычисления производных от f(x,y), а требуют только вычисления
    самой функции.

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

1.1. Метод Эйлера.

Этот метод, один из самых старых и широко известных, описывается формулой:

        y[m+1] = y[m] + h*f(x[m],y[m]).                 (3)


Найденное по формуле (3) решение согласуется с разложением в ряд Тейлора
вплоть до членов порядка h, т.е. данный метод является методом Рунге-Кутта
первого порядка.
Этот метод имеет довольно большую ошибку приближения, кроме того, он очень
часто оказывается неустойчивым - изначально малая ошибка (происходящая от
приближения, округления или исходных данных) увеличивается с ростом x.
Для вычисления y[m+1] метод Эйлера использует наклон касательной только в
точке x[m],y[m]. Этот метод можно усовершенствовать множеством различных
способов. Рассмотрим два из них.

1.2. Исправленный метод Эйлера.

В исправленном методе Эйлера мы находим средний tg угла наклона касательной
для двух точек: [I]x[m],y[m] и x[m+1],y[m]+h*y'[m][/I]. Соотношения, описывающие
данный метод, имеют вид:

        y[m+1] = y[m] + h*F(x[m],y[m],h)                (4)
        F(x[m],y[m],h) = ( y'[m] + f(x[m]+h, y[m]+h*y'[m]) )/2      (5)
        y'[m] = f(x[m],y[m])                        (6)


Исправленный метод Эйлера согласуется с разложением в ряд Тейлора вплоть до
членов степени h^2, являясь, т.о. методом Рунге-Кутта второго порядка.

1.3. Модифицированный метод Эйлера.

В данном методе мы находим tg угла наклона касательной в точке:

        x = x[m] + h/2;     y = y[m] + (h/2)*y'[m]


Соотношения, описывающие модифицированный метод Эйлера имеют вид:

        y[m+1] = y[m] + h*F(x[m],y[m],h)                (7)
        F(x[m],y[m],h) = f( x[m]+h/2, y[m]+(h/2)*y'[m] )        (8)
        y'[m] = f(x[m],y[m])                        (9)


Модифицированный метод Эйлера также согласуется с разложением в ряд Тейлора
вплоть до членов степени h^2, и также является методом Рунге-Кутта второго
порядка.

1.4. Метод Рунге-Кутта четвертого порядка.

Данный метод является одним из самых употребительных методов интегрирования
дифференциальных уравнений. Этот метод применяется настолько широко, что в
литературе его просто называют "методом Рунге-Кутта" без всяких указаний на его
тип или порядок. Этот классический метод описывается системой следующих пяти
соотношений:

        y[m+1] = y[m] + h*(k[1] + 2*k[2] + 2*k[3] + k[4])/6     (10)
        k[1] = f(x[m], y[m])                        (11)
        k[2] = f(x[m]+h/2, y[m]+h*k[1]/2)               (12)
        k[3] = f(x[m]+h/2, y[m]+h*k[2]/2)               (13)
        k[4] = f(x[m]+h, y[m]+h*k[3])                   (14)


Ошибка приближения для этого метода равна e[t]=k*h^5. Заметим, что при
использовании этого метода функцию необходимо вычислять четыре раза.

Один из серьезных недостатков методов Рунге-Кутта состоит в отсутствии
простых способов оценки их ошибки.

2. Методы прогноза и коррекции.

Отличительной чертой методов Рунге-Кутта является то, что при вычислении
следующей точки (x[m+1],y[m+1]) используется информация только об одной
предыдущей точке (x[m],y[m]), но не нескольких. Кроме того, для методов
Рунге-Кутта отсутствует достаточно простые способы оценки ошибки, что приводит
к необходимости рассмотрения некоторых дополнительных методов решения ДУ.
Отличительное свойство этих методов состоит в том, что с их помощью нельзя
начать решение уравнения т.к. в них необходимо использовать информацию о
предыдущих точках решения. Чтобы начать решение, имея только одну точку,
определяемую начальными условиями, или для того, чтобы изменить шаг (h),
необходим метод типа Рунге-Кутта. Поэтому приходится использовать разумное
сочетание этих двух методов.
Методы, которые мы рассмотрим, известны под общим названием методов прогно-
за и корректировки. Как ясно из названия вначале "предсказывается" значение
y[m+1], а затем используется тот или иной метод его "корректировки".
Среди множества возможных формул прогноза и коррекции выберем по одному
примеру, применимому ко многим практическим задачам.

2.1. Метод Адамса-Бошфора.

Для прогноза используем формулу второго порядка

        y^(0)_[m+1] = y[m-1] + 2*h*f(x[m],y[m])             (15)


где (0) - означает исходное приближение y[m+1], т.е. предсказанное значение.
Непосредственно из написанной формулы следует, что с ее помощью нельзя вычис-
лить y[1], т.к. для его вычисления потребовалась бы точка, расположенная перед
начальной точкой y[0]. Чтобы начать решение с помошью метода прогноза и коррек-
ции, для нахождения y[1] необходимо использовать метод типа Рунге-Кутта.
Для коррекции возьмем формулу, похожую на исправленный метод Эйлера (4)-(6):

        y^(i)_[m+1] = y[m] + h(f(x[m],y[m])+f(x[m+1],y^(i-1)_[m+1]))/2  (16)


для i=1,2,3, ...
Итерационный процесс прекращается, когда
        | y^(i+1)_[m+1] - y^(i)_[m+1] | < Eps               (17)


для некоторого Eps>0.

    {
      Writeln('Labinskiy Nikolay aka e-moe (c) 2005');
      Writeln('See you later on forum.sources.ru ;)');
    }
     
    {$E+,N+,X+}
    const
      eps = 0.00001;
      x0  = 3;
      y0  = 3;
     
    type
      TFunc = function(x,y: extended): extended;
     
    function func(x,y: extended): extended; far;
    begin
      func := y*cos(x) - 2*sin(2*x);
    end;
     
    function Euler(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
    { где x0_ и y0_  - начальное условие
          x_end - точка, в которой необходимо вычислить результат
          n - количество шагов для вычисления результата }
    var
      i   : word;      { счетчик цикла }
      x,h : extended;  { текущая точка и длина шага }
      res : extended;  { переменная для накопления конечного результата функции}
    begin
      h:= (x_end - x0_)/n; { Находим длину шага }
      res:= y0_; { устанавливаем начальные значения}
      x:=x0_;
      for i:=1 to n do
        begin { вычисляем результат по методу Эйлера }
          res:=res+h*f(x,res);
          x:=x+h; { переходим к следующей точке }
        end;
      Euler:=res;  { присваиваем конечный результат функции }
    end;
     
    function Euler2(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
    { где x0_ и y0_  - начальное условие
          x_end - точка, в которой необходимо вычислить результат
          n - количество шагов для вычисления результата }
    var
      i   : word;      { счетчик цикла }
      x,h : extended;  { текущая точка и длина шага }
      res : extended;  { переменная для накопления конечного результата функции}
    begin
      h:= (x_end - x0_)/n; { Находим длину шага }
      res:= y0_; { устанавливаем начальные значения}
      x:=x0_;
      for i:=1 to n do
        begin { вычисляем результат по исправленному методу Эйлера }
          res:=res+h*(f(x,res)+f(x+h,res+h*f(x,res)))/2;
          x:=x+h; { переходим к следующей точке }
        end;
      Euler2:=res; { присваиваем конечный результат функции }
    end;
     
    function Euler3(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
    { где x0_ и y0_  - начальное условие
          x_end - точка, в которой необходимо вычислить результат
          n - количество шагов для вычисления результата }
    var
      i   : word;      { счетчик цикла }
      x,h : extended;  { текущая точка и длина шага }
      res : extended;  { переменная для накопления конечного результата функции}
    begin
      h:= (x_end - x0_)/n; { Находим длину шага }
      res:= y0_; { устанавливаем начальные значения}
      x:=x0_;
      for i:=1 to n do
        begin { вычисляем результат по модифицированному методу Эйлера }
          res:=res+h*f(x+h/2,res+(h/2)*f(x,res));
          x:=x+h; { переходим к следующей точке }
        end;
      Euler3:=res; { присваиваем конечный результат функции }
    end;
     
    function RungeKutt(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
    { где x0_ и y0_  - начальное условие
          x_end - точка, в которой необходимо вычислить результат
          n - количество шагов для вычисления результата }
    var
      i   : word;      { счетчик цикла }
      x,h : extended;  { текущая точка и длина шага }
      res : extended;  { переменная для накопления конечного результата функции }
      k1,k2,k3,k4: extended; { вспомогательные переменные вычисления результата }
    begin
      h:= (x_end - x0_)/n; { Находим длину шага }
      res:= y0_; { устанавливаем начальные значения}
      x:=x0_;
      for i:=1 to n do
        begin { вычисляем результат по методу Рунге-Кутта 4го порядка }
          k1:=f(x,res);
          k2:=f(x+h/2,res+h*k1/2);
          k3:=f(x+h/2,res+h*k2/2);
          k4:=f(x+h,res+h*k3);
          res:=res+h*(k1+2*k2+2*k3+k4)/6;
          x:=x+h; { переходим к следующей точке }
        end;
      RungeKutt:=res; { присваиваем конечный результат функции }
    end;
     
    function AdmsBoshf(f: TFunc; x0_,x_end,y0_: extended; n: word): extended;
    { где x0_ и y0_  - начальное условие
          x_end - точка, в которой необходимо вычислить результат
          n - количество шагов для вычисления результата }
    var
      i     : word;         { счетчик цикла }
      x,h   : extended;     { текущая точка и длина шага }
      y1,y2 : extended;     { переменные для вычисления следующей точки }
      tmp   : extended;     { временная переменная для уточнения результата }
    begin
      h:= (x_end - x0_)/n; { Находим длину шага }
      x:=x0_; { устанавливаем текущую точку }
      y1:=y0_+h*f(x+h/2,y0_+(h/2)*f(x,y0_)); { находим начальное значение
                                               модифицированным методом Эйлера }
      repeat { уточняем(корректируем) начальное значение }
        tmp:=y1;
        y1:=y0_+h*(f(x,y0_)+f(x+h,y1))/2;
      until abs(y1-tmp) < eps;
     
      x:=x+h; { переходим к следующей точке }
      y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение }
     
      repeat { уточнаяем наш прогноз }
        tmp:=y2;
        y2:=y1+h*(f(x,y1)+f(x+h,y2))/2;
      until abs(y2-tmp) < eps;
     
      for i:=3 to n do { начинаем счет с 3 т.к. две точки уже найдены }
        begin
          x0_:=x; { переходим к следующей точке }
          x:=x+h;
          y0_:=y1;
          y1:=y2;
     
          y2:=y0_+2*h*f(x,y1); { прогнозируенм следующее значение }
     
          repeat { уточнаяем наш прогноз }
            tmp:=y2;
            y2:=y1+h*(f(x,y1)+f(x+h,y2))/2;
          until abs(y2-tmp) < eps;
        end;
     
      AdmsBoshf:=y2; { присваиваем конечный результат функции }
    end;
     
    begin
      writeln('Численное решение дифференциальных уравнений:');
      writeln(#10,'   y'' = y*cos(x) - 2*sin(2*x);   y(0)=3;   x_end=',(5*x0+3.5):5:5);
      writeln(#10,'Метод Эйлера:');
      writeln('n=5 000, result: ',Euler(func,x0,5*x0+3.5,y0,5000):5:5);
      writeln('n=10 000, result: ',Euler(func,x0,5*x0+3.5,y0,10000):5:5);
      writeln('n=25 000, result: ',Euler(func,x0,5*x0+3.5,y0,25000):5:5);
      writeln(#10,'Исправленный метод Эйлера:');
      writeln('n=50, result: ',Euler2(func,x0,5*x0+3.5,y0,50):5:5);
      writeln('n=100, result: ',Euler2(func,x0,5*x0+3.5,y0,100):5:5);
      writeln('n=250, result: ',Euler2(func,x0,5*x0+3.5,y0,250):5:5);
      writeln(#10,'Модифицированный метод Эйлера:');
      writeln('n=50, result: ',Euler3(func,x0,5*x0+3.5,y0,50):5:5);
      writeln('n=100, result: ',Euler3(func,x0,5*x0+3.5,y0,100):5:5);
      writeln('n=250, result: ',Euler3(func,x0,5*x0+3.5,y0,250):5:5);
     
      writeln(#10,'Метод Рунге-Кутта:');
      writeln('n=50, result: ',RungeKutt(func,x0,5*x0+3.5,y0,50):5:5);
      writeln('n=100, result: ',RungeKutt(func,x0,5*x0+3.5,y0,100):5:5);
      writeln('n=250, result: ',RungeKutt(func,x0,5*x0+3.5,y0,250):5:5);
     
      write(#10,'Press Enter to continue.');
      readln;
     
      writeln(#10,'Метод Адамса-Бошфора:');
      writeln('n=50, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,50):5:5);
      writeln('n=100, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,100):5:5);
      writeln('n=250, result: ',AdmsBoshf(func,x0,5*x0+3.5,y0,250):5:5);
    end.


Автор статьи: Labinskiy Nikolay aka e-moe © 2005

Это сообщение было перенесено сюда или объединено из темы "Численное решение ДУ"]]> Romtek Pascal: Математика https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=670695 Mon, 04 Apr 2005 20:38:08 +0000 Численные методы https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=670695 Численное решение систем линейных алгебраических уравнений

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

1. Метод исключений (метод Гаусса).

Данный метод является наиболее известным и широко применяемым.

Рассмотрим систему из n уравнений с n неизвестными. Обозначим неизвестные через x[1], x[2], ... , x[n] и запишем систему в следующем виде:

         a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
         a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
         ....................................................
         a[n,1]*x[1] + a[n,2]*x[2] + ... + a[n,n]*x[n] = b[n]


Предполагается, что в силу расположения уравнений, a[i,i] <> 0. Если это не так, то, меняя уравнения местами, добиваемся выполнения этого условия. Введем n-1 множителей:

         m[i] = a[i,1]/a[1,1]   , где i = 2, 3, ... , n


и вычтем из каждого i-го уравнения первое, умноженное на m, обозначая

         a'[i,j] = a[i,j] - m[i]*a[1,j]
         b'[i] = - m[i]*b[i]
         i = 2, 3, ... , n;  j = 1, 2, ... , n.


Для всех уравнений, начиная со второго, получим:

         a[i,1] = 0   , где i = 2, 3, ... , n.


Преобразованная система уравнений запишется в следующем виде:

         a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
            0   +     a'[2,2]*x[2] + ... + a'[2,n]*x[n] = b'[2]
         ....................................................
            0   +     a'[n,2]*x[2] + ... + a'[n,n]*x[n] = b'[n]


Продолжая таким же образом, исключаем x[2] из последних n-2 уравнений, затем x[3] из последних n-3 и т.д. На некотором k-м этапе мы исключаем x[k] с помощью множителей:

         m^(k-1)_[i] = a^(k-1)_[i,k]/a^(k-1)_[k,k]   , i = k+1, ... , n,


не забывая, что a(k-1)[k,k] <> 0.
Примечание, ^(k-1) обозначает на k-1й итерации
_[i,k] - индекс эл-та.
Тогда

         a^(k)_[i,j] = a^(k-1)_[i,j] - m^(k-1)_[i]*a^(k-1)_[k,j]
         b^(k)_[i] = - m^(k-1)_[i]*b^(k-1)_[k]
         i = k+1, k+2, ... , n;  j = k, k+1, ... , n.


Окончательно, треугольная система уравнений записывается:

         a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
                       a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
         ....................................................
                                           a[n,n]*x[n] = b[n]


Обратная подстановка для нахождения значений неизвестных задается формулами:

         x[n] = b^(n-1)_[n] / a^(n-1)_[n,n]
         x[n-1] = (b^(n-2)_[n-1] - a^(n-2)_[n-1,n]*x[n]) / a^(n-2)_[n-1,n-1]
         ...................................................................
         x[j] = ( b^(j-1)_[j] - a^(j-1)_[j,n]*x[n] - ... -
                - a^(j-1)_[j,j+1]*x[j+1] ) / a^(j-1)_[j,j]
     
         j = n-2, n-3, ... , 1.


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

         | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] |


Поэтому перед началом процесса исключения очередного неизвестного необходимо переставить уравнения так, чтобы максимальный по модулю коэффициент при x[k] попал на главную диагональ (данный способ решения часто называется методом главного элемента).

2. Итерационный метод Гаусса-Зейделя.

Данный метод отличается простотой и легкостью программирования, малой ошибкой округления, но сходится при выполнении некоторых условий.

Рассмотрим все ту же систему из n уравнений с n неизвестными. По прежнему полагаем, что a[i,i] <>0 для всех i. Тогда k-е приближение к решению будет задаваться формулой:

      x^(k)_[i] = ( b[i] - a[i,1]*x^(k)_[1] - ... - a[i,j-1]*x^(k)_[i-1] -
                  - a[i,j+1]*x^(k-1)_[i+1] - a[i,n]*x^(k-1)_[n] ) / a[i,i]
      i = 1, 2, ... , n.


Как известно, любой итерационный процесс, как и настоящий мужик, должен имееть свой конец :)

Поэтому вычисления нужно продолжать до тех пор, все x^(k)_[i] не станут достаточно близки к x(k-1)[i], т.е. пока для всех i не будет выполнятся неравенство

         max| x^(k)_[i] - x^(k-1)_[i] | < Eps,


где Eps некоторое положительное число, задающее точность вычислений.

Итерационный метод Гаусса-Зейделя сходится, если удовлетворяется достаточное условие сходимости: для всех i модуль диагонального коэффициента должен быть не меньше суммы модулей остальных коэффициентов данной строки

         | a[i,i] | >= | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | +
                       + | a[i,j+1] | + ... + | a[i,n] | ,


а хотя бы для одной строки это неравенство должно выполнятся строго

         | a[i,i] | > | a[i,1] | + | a[i,2] | + ... + | a[i,j-1] | +
                      + | a[i,j+1] | + ... + | a[i,n] |.


3. Сравнение методов.

Метод исключений имеет то преимущество, что он конечен и теоретически, с его помощью можно решить любую невырожденную систему. Итерационный метод сходится только для специальных систем уравнений, однако когда итерационные методы сходятся, они, обычно, предпочтительнее:
  • Время вычислений пропорционально n2 на итерацию, в то время как для метода исключений время пропорционально n3. Если для решения системы требуется меньше чем n итераций, то общие затраты машинного времени будут меньше.
  • Как правило, ошибки округления в итерационном методе меньше. Часто, это соображение может оказаться достаточно важным, что оправдывает дополнительные затраты машинного времени. Многие систему, возникающие на практике, имеют среди коэффициентов большое кол-во нулей. В этом случае итерационные методы, если они сходятся, в высшей степени предпочтительны т.к. в методе исключений получается треугольная система, которая в качестве коэффициентов нулей обычно не содержит.

А вот и исходник по теме:

    {
      Writeln('Labinskiy Nikolay aka e-moe (c) 2005');
      Writeln('See you later on forum.sources.ru ;)');
    }
     
    {$N+,G+,X+}
    const
      Comments: boolean = false; { нужно ли печатать пошаговые рез-ты расчета }
      Eps = 0.00001;        { Точность вычислений }
      n_max = 4;            { Кол-во ур-й и неизвестных }
      const_AnB: array[1..n_max,1..n_max+1] of double =
     (( -8,   2,  17,  -5, -119.97),
     { -8*X[1] + 2*X[2] + 17*X[3] - 5*X[4] = -119.97  и т.д. ... }
      (  4, -22,   6,   5,   55.73),
      ( 15,   3,  -5,  -5,   18.77),
      ( -4,  -4,   5,  14,  -79.42));
    { Описание типов матрицы уравнеиний и массива найденных Х }
    Type TMatr = array[1..n_max,1..n_max+1] of double;
    Type TXMarr = array[1..n_max] of double;
     
    procedure PrintX(const X: TXMarr; n: byte; cnt: word);
    { Печатает найденные Х }
    var
      k,       { Позиция заданного корня в массиве }
      i: byte; { Номер текущего корня для вывода }
    begin
      if n = 0 then
        exit;
      k:=low(X);
      if cnt = 0 then
        begin
          write(' N       x1');
          for i:=2 to n do
            write('         x',i);
          writeln;
        end;
      write(cnt:2);
      for i:=1 to n do
        begin
          Write('   ',X[k]:3:5);
          inc(k);
        end;
      writeln;
    {
      При решении систем, с большим кол-вом неизвестных, возможно,
      придется переделать эту процедуру для корректного вывода
      всех решений.
    }
    end; { PrintX }
     
    procedure PrintMatr(var AnB: TMatr; n: byte);
    { Печатает заданную матрицу }
    var
      i,j: byte;
    begin
      for i:=1 to n do
        begin
          for j:=1 to n+1 do
            write(AnB[i,j]:3:5,'   ');
          writeln;
        end;
      writeln;
    {
      При решении систем, с большим кол-вом неизвестных, возможно,
      придется переделать эту процедуру для корректного вывода
      всех решений.
    }
    end; { PrintMatr }
     
    function Normalize(var AnB: TMatr; n: byte): byte;
    {
     Располагает на главной диагонали наибольшие элементы в столбцах
     Если один из них равен нулю,то систему решить не получится ...
     Возвращает 0 если систему заданными методами решить не получится,
     1 - с-му можно решить только методом Гаусса
     2 - с-му можно решать любым методом
    }
    var
      max: double; { Переменная для поиска макс. эл-та }
      imax: byte;  { Переменная для поиска макс. эл-та }
      tmp: double; { Переменная для обмена элементов }
      i,j: byte;   { Счетчики циклов }
    begin
      Normalize:=2;
      for j:=1 to n do
        begin
          max:=Abs(AnB[1,j]);
          imax:=1;
          for i:=2 to n do
            if Abs(AnB[i,j]) > max then
              begin
                max:=Abs(AnB[i,j]);
                imax:=i;
              end;
          if max < Eps then
            begin
              Normalize:=0;
              Writeln('Error: a[i,i]=0!');
              exit;
            end
          else
            if j <> imax then
              for i:=1 to n+1 do
                begin
                  tmp:=AnB[j,i];
                  AnB[j,i]:=AnB[imax,i];
                  AnB[imax,i]:=tmp;
                end;
        end;
      for i:=1 to n do
        begin
          tmp:=0;
          for j:=1 to n do
            if i <> j then
              tmp:=tmp+Abs(AnB[i,j]);
          if Abs(AnB[i,i]) < tmp then
            Normalize:=1;
        end;
    end; { Normalize }
     
    procedure Metod1(const Matr; n: byte);
    var
      k    : byte;
      M_   : TMatr absolute Matr;
      AnB  : TMatr;
      X,M  : TXMarr;
    function Calc(k: byte): boolean;
    var
      i,j  : byte;
    begin
      Calc:= true;
      for i:=k+1 to n do
        begin
          M[i]:= AnB[i,k]/AnB[k,k];
          if M[i] > 1 then
            begin
              Calc:= false;
              exit;
            end;
        end;
     
      for i:=k+1 to n do
        for j:=k to n+1 do
          AnB[i,j]:=AnB[i,j]-M[i]*AnB[k,j];
     
      if Comments then
        PrintMatr(AnB,n);
    end; { Calc }
    function Summ(j: byte): double;
    var
      i : byte;
      res : double;
    begin
      res:=AnB[j,n+1];
      for i:=n downto j+1 do
        res:=res - AnB[j,i]*X[i];
      Summ:=res;
    end; { Summ }
    begin
      Writeln('Метод Гаусса:');
      AnB:=M_;
      if Normalize(AnB,n) = 0 then
        begin
          writeln('Error: Невозможно решить систему уравнений!');
          exit;
        end;
     
      for k:=1 to n-1 do
        if not Calc(k) then
          begin
            Writeln('Error: Невыполняется условие | a^(k-1)_[k,k] | >= | a^(k-1)_[i,k] | !');
            exit;
          end;
     
      X[n]:=AnB[n,n+1]/AnB[n,n];
      for k:=n-1 downto 1 do
        X[k]:=Summ(k)/AnB[k,k];
     
      PrintX(X,n,0);
      Writeln;
    end; { Metod1 }
     
    procedure Metod2(const Matr; n: byte);
    var
      M_   : TMatr absolute Matr;
      AnB  : TMatr;
      X    : TXMarr;
      i,j: byte;
      tmp: double;
      delta: double;
      cnt: word;
    function Summ(i: byte): double;
    var
      j: byte;
      res: double;
    begin
      res:=AnB[i,n+1];
      for j:=1 to n do
        if i <> j then
          res:=res-AnB[i,j]*X[j];
      Summ:=res;
    end; { Summ }
    begin
      writeln('Метод Гаусса-Зейделя:');
      AnB:=M_;
      For i:=1 to n do
        X[i]:=0;
      if Normalize(AnB,n) <> 2 then
        begin
          writeln('Error: !');
          exit;
        end;
      delta:=1;
      cnt:=0;
      while delta > Eps do
        begin
          if Comments then
            PrintX(X,n,cnt);
          delta:=0;
          for i:=1 to n do
            begin
              tmp:=Summ(i)/AnB[i,i];
              if delta < Abs(tmp-X[i]) then
                delta:=Abs(tmp-X[i]);
              X[i]:=tmp;
            end;
          inc(cnt);
        end;
     
      PrintX(X,n,cnt);
     
      writeln('Заданная точность вычислений достигнута на ',cnt,' шаге.');
      Writeln;
    end; { Metod2 }
     
    begin
      Writeln(#10,'Численное решение систем линейных алгебраических уравнений:',#10);
      Metod1(const_AnB,n_max);
      Writeln('Press Enter to continue');
      readln;
      Metod2(const_AnB,n_max);
      Writeln('Press Enter if you wanna DOS ;)');
      readln;
    end.


Это сообщение было перенесено сюда или объединено из темы "Решение уравнений"]]> e-moe Pascal: Математика https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=665661 Thu, 31 Mar 2005 11:59:35 +0000 Численные методы https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=665661 Численное решение уравнений итерационными методами.

Автор:
  DonNTU. e-moe aka Labinskiy Nikolay © 2005


Данная статья посвящена нахождению корней уравнения
            F(x)=0,    (1)

где функция F(x) может быть алгебраической либо трансцендентной и должна удовлетворять условию дифференцируемости.
Как правило, численное решение уравнений состоит из двух этапов:
нахождение приближенного значения корня и его уточнение до заданной точности. Начальное приближение часто
известно из физических соображений либо находится специальными методами, например, графически. Подробно будет
рассмотрен второй этап решения уравнений: нахождение корня с заданной точностью различными итерационными методами.
  • Метод последовательных приближений.
      function Metod1(f: TFunc; x0: double): double;

    В данном методе для удобства вычислений переходят от исходного уравнения (1) к уравнению
              x=f(x).                                         (2)

    Данный переход можно осуществить множеством способов, например, прибавить к обеим частям (1) x.
    Суть метода последовательных приближений состоит в том, что начальное приближение x[0] подставляется в правую часть
    формулы (2) и вычисляется значение x[1]. Затем полученное x[1] снова подставляется в правую часть формулы (2)
    и вычисляется x[2], потом x[3] и т.д. Рабочая формула метода последовательных приближений имеет вид
              x[n]=f(x[n-1]).                                 (3)

    Вычисления продолжаются до тех пор, пока не будет достигнута заданная точность Eps, т.е
              | x[n]-x[n-1] | < Eps.                          (4)

    Основной проблемой при работе с итерационными методами является проблема сходимости. Достаточным условием сходимости метода последовательных приближений является выполнение условия
              | f'(x[n]) | < 1                                (5)
    для всех значений x[n].
  • Усовершенствованный метод последовательных приближений.
      function Metod2(f: TFunc; x0: double): double;

    Формула данного итерационного метода имеет вид
              x[n+1] = x[n] + A*( f(x[n]) - x[n] ),           (6)
    где A определяется по формулам
              A = 1/(1 - f'(s))                               (7)
              f'(s) = ( f(x[n]) - x[n] )/( x[n] - x[n-1] )    (8)
    при этом на первом шаге x[1] определяется простым методом последовательных приближений.
  • Метод Ньютона-Рафсона, Бриге-Виетта.
      function Metod3(f: TFunc; x0: double): double;

    Небольшая дальнейшая модификация усовершенствованного метода последовательных приближений приводит к одному из
    наиболее известных численных методов решения уравнений - методу Ньютона-Рафсона. Формула метода для f(x), подчиняющегося соотношению (2), имеет вид
              x[n+1] = (f(x[n])-x[n]*f'(x[n]))/(1-f'(x[n])),  (9)
    при этом сходимость метода обеспечивается, если:
    1. x[0] выбрано достаточно близко к решению x = f(x);
    2. Вторая производная f''(x) не становится слишком большой;
    3. Производная f'(x) не слишком близка к 1.
      Формула Ньютона-Рафсона для F(x), подчиняющегося соотношению (1), имеет вид
                x[n+1] = x[n] - F(x[n])/F'(x[n]),               (10)
при этом условия сходимости принимают вид:
  1. x[0] выбрано достаточно близко к решению x = F(x);
  2. Вторая производная f''(x) не становится слишком большой;
  3. Производная f'(x) не слишком близка к 0.

Применим метод Ньютона-Рафсона согласно формуле (10), при этом вычисление F(x) будем осуществлять по правилу Горнера
с использованием рекуррентных формул:
            b[m] = a[m]
            b[j] = a[j] + x[n]*b[j+1], j=m-1, ..., 0        (11)

Таким образом находим F(x) = b[0].
F'(x) представляет собой многочлен степени m-1. Воспользовавшись для его вычисления теми же рекуррентными формулами, имеем:
            c[m] = b[m]
            c[j] = b[j] + x[n]*c[j+1], j=m-1, ..., 1        (12)

и соответственно F'(x) = c[1].
Подставляя найденные значения F(x) и F'(x) в формулу (10) для метода Ньютона-Рафсона, получаем
            x[n+1] = x[n] - b[0]/c[1],                      (13)

где b[o] и c[1] вычислены по формулам (11) и (12).
Тихо и незаметно мы пришли к методу Бриге-Виетта ;)
    {
      DonNTU. e-moe aka Labinskiy Nikolay (c) 2005
      visit www.sources.ru
    }
     
    {$N+,G+,X+}
    const
      Eps = 0.00001;        { Точность вычислений }
      anMax = 3;            { Степень полинома }
      { Коэффициенты полинома F(x) = a[0]+a[1]*x+a[2]*x^2+ ... +a[n]*x^n }
      An: array[0..anMax] of double = (-5.372,1.2493,0.559,-0.13);
                                      { a[0]   a[1]  a[2]  a[3] }
    type
      { Описание функции для передачи ее в кач-ве параметра функциям }
      TFunc = function(x: double; px: byte; d: shortint): double;
    var
      x0,x1,x2: double; { Переменные для вычисления корней }
     
    function Func(x: double; px: byte; d: shortint): double; far;
    (*
      Функция для вычисления полинома в точке 'x'
      x  - в этой точке будет происходить расчет
      px - если 1, то результат будет = F(x)+x
      d  - если 0, то F(x), 1 - F'(x), -1 - F(x)/F'(x)
    *)
    var
      bn,                   { Это будет F(x) }
      cn: double;           { Это будет F'(x) }
      j: byte;              { Счетчик цикла ;) }
    begin
      if px=1 then          { Если px=1, то добавляем к рез-ту x }
        An[1]:=An[1]+1;
    { Далее идет расчет по правилу Горнера  }
      bn:=An[anMax];
      cn:=bn;
      for j:= anMax-1 downto 1 do
        begin
          bn:=An[j]+x*bn;
          cn:=bn+x*cn;
        end;
      bn:=An[0]+x*bn;
    { В зависимости от d возвращаем результат }
      case d of
        -1: Func:=bn/cn;
         0: Func:=bn;
         1: Func:=cn;
      else
        Func:=bn;
      end;
    { Забираем обратно свой x }
      if px=1 then
        An[1]:=An[1]-1;
    end;
     
    function Sign(x: double): shortint;
    { Функция определения знака }
    begin
      if x < - Eps then
        Sign:=-1
      else
        if x > Eps then
          Sign:=1
        else
          Sign:=0;
    end;
     
    function Get_X0(f: TFunc): double;
    (*
      Поиск начального приближения:
      Идет начальный просмотр функции и определение двух точек,
      между которыми функция меняят знак (проходит через 0).
      Где-то между этими точками лежит корень.
    *)
    const
      LookFrom = -10.0;     { Поиск начинается с этой точки }
      step = 1;             { Шаг, с которым производится поиск }
    begin
      x0:=LookFrom;
      x1:=x0+eps;
      while Sign(f(x1,0,0)) = Sign(f(x0,0,0)) do
        begin
          x0:=x1;
          x1:=x1+step;
        end;
      Get_X0:= (x1+x0)/2;
    end;
     
    function Metod1(f: TFunc; x0: double): double;
    var
      cnt: word;
    function TestFunc(x: double): boolean;
    begin
      if abs(f(x,1,1)) >= 1 then
        begin
          writeln('  Ошибка: Невозможно расчитать результат!');
          writeln('  Модуль производной в точке ',x:0:6,' больше 1',#10);
          TestFunc:= false;
        end
      else
        TestFunc:= true;
    end;
    begin
      Writeln('Метод последовательных приближений:');
      if TestFunc(x0) then
        x1:=f(x0,1,0)
      else
        exit;
    {  x2:=f(x1,0)+x1;}
      cnt:=1;
      while abs(x1-x0) > eps do
        begin
          x0:=x1;
          if TestFunc(x0) then
            x1:=f(x0,1,0)
          else
            exit;
          inc(cnt)
        end;
      Metod1:=x1;
      Writeln('x=', x1:0:6);
      Writeln('Метод сошелся на ',cnt,' шаге.',#10);
    end;
     
    function Metod2(f: TFunc; x0: double): double;
    var
      alph: double;
      cnt: word;
    begin
      Writeln('Усовершенствованный метод последовательных приближений:');
      x1:=f(x0,1,0);
      alph:=1/(1-f(x1,0,0)/(x1-x0));
      x2:=x1+alph*f(x1,0,0);
      cnt:=2;
      while abs(x2-x1) > eps do
        begin
          x0:=x1;
          x1:=x2;
          alph:=1/(1-f(x1,0,0)/(x1-x0));
          x2:=x1+alph*f(x1,0,0);
          inc(cnt)
        end;
      Metod2:=x2;
      Writeln('x=', x2:0:6);
      Writeln('Метод сошелся на ',cnt,' шаге.',#10);
    end;
     
    function Metod3(f: TFunc; x0: double): double;
    var
      cnt: word;
    begin
      Writeln('Метод Бирге-Виетта:');
      x1:=x0-f(x0,0,-1);
      cnt:=1;
      while abs(x1-x0) > Eps do
        begin
          x0:=x1;
          x1:=x0-f(x0,0,-1);
          inc(cnt);
        end;
      Metod3:=x1;
      Writeln('x=', x1:0:6);
      Writeln('Метод сошелся на ',cnt,' шаге.',#10);
    end;
     
    begin
      x0:=Get_X0(Func);
      Writeln;
      Writeln('F(x)= -5.372 + 1.2493*x + 0.559*x^2 - 0.13*x^3');
      Writeln('Eps = ',eps:0:5,'  x0= ',x0:0:6,#10);
     
      Metod1(Func,x0);
      Metod2(Func,x0);
      Metod3(Func,x0);
     
      Write('Press Enter to exit');
      Readln;
    end.


Это сообщение было перенесено сюда или объединено из темы "Решение уравнений"]]> Romtek Pascal: Математика https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=653753 Tue, 22 Mar 2005 09:54:08 +0000 Численные методы https://forum.sources.ru/index.php?showtopic=100241&view=findpost&p=653753 Решение уравнений

    {Автор: Ozzя}
     
    var
      a,b,x,e:real;
     
    { Нелинейное уравнение (x-2)*lg(x+11)-1 }
    function f(x:real):real;
    begin
      f:=(x-2)*lg(x+11)-1;
    end;
     
    { Функция, определяющая знак числа }
    function sgn(x:real):integer;
    begin
      sgn:=0;
      if x<0.0 then
        sgn:=-1;
      if x>0.0 then
        sgn:=1;
    end;
     
    { Процедура, реализующая метод половинного деления }
    procedure dich (var a,b,x,e:real);
    var
      i:integer;
      r:real;
    begin
      i:=sgn(f(a));
      { Повторяем до тех пор, пока интервал [a,b] не станет меньше   }
      { заданной погрешности е                                       }
      while b-a>e do
        begin
          { Определяем середину отрезка                              }
          x:=(a+b)/2;
          { Далее делаем выбор, какую из частей отрезка взять        }
          { для дальнейшего уточнения корня                          }
          { Если левая часть уравнения f(x) есть непрерывная функция }
          { аргумента х, то корень будет находиться в той половине   }
          { отрезка, на концах которой f(x) имеет разные знаки.      }
          r:=f(x);
          if sgn(r)=i then
            a:=x
          else
            b:=x
        end
    end;
     
    begin
      write('Введите значения интервала A,B? ');
      readln(a,b);
      write('Введите требуемое значение точности E? ');
      readln(e);
      dich(a,b,x,e);
      writeln('X= ',x);
    end.


Это сообщение было перенесено сюда или объединено из темы "Решение уравнений"

Это сообщение было перенесено сюда или объединено из темы "Заготовка для "Численных методов""]]> Romtek Pascal: Математика