На главную
ПРАВИЛА FAQ Помощь Участники Календарь Избранное DigiMania RSS
msm.ru
! Правила ЧаВо (FAQ) разделов Паскаля
В этом разделе разрешено создавать только темы, в которых описано РЕШЕНИЕ какой-либо общей проблемы, или описание какого-либо аспекта языка Паскаль.
Обсуждение уже созданных тем разрешено, но только конструктивное, например указание на ошибку или уточнение имеющегося текста.

Также читать Требования к оформлению статей
Модераторы: Romtek, volvo877
  
    > Численные методы, Уравнения, СЛАУ, пределы, ряды
      Решение уравнений

      ExpandedWrap disabled
        {Автор: 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.


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

      Это сообщение было перенесено сюда или объединено из темы "Заготовка для "Численных методов""
        Численное решение уравнений итерационными методами.

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


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

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

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

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

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

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

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

          Небольшая дальнейшая модификация усовершенствованного метода последовательных приближений приводит к одному из
          наиболее известных численных методов решения уравнений - методу Ньютона-Рафсона. Формула метода для f(x), подчиняющегося соотношению (2), имеет вид
          ExpandedWrap disabled
                    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), имеет вид
            ExpandedWrap disabled
                      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) будем осуществлять по правилу Горнера
        с использованием рекуррентных формул:
        ExpandedWrap disabled
                  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. Воспользовавшись для его вычисления теми же рекуррентными формулами, имеем:
        ExpandedWrap disabled
                  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) для метода Ньютона-Рафсона, получаем
        ExpandedWrap disabled
                  x[n+1] = x[n] - b[0]/c[1],                      (13)

        где b[o] и c[1] вычислены по формулам (11) и (12).
        Тихо и незаметно мы пришли к методу Бриге-Виетта ;)
        ExpandedWrap disabled
          {
            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 -
          Численное решение систем линейных алгебраических уравнений


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

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

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

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

          ExpandedWrap disabled
                 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 множителей:

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


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

          ExpandedWrap disabled
                 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.


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

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


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

          ExpandedWrap disabled
                 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] с помощью множителей:

          ExpandedWrap disabled
                 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] - индекс эл-та.
          Тогда

          ExpandedWrap disabled
                 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.


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

          ExpandedWrap disabled
                 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]


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

          ExpandedWrap disabled
                 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.


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

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


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

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

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

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

          ExpandedWrap disabled
              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 не будет выполнятся неравенство

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


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

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

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


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

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


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

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

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

          ExpandedWrap disabled
            {
              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.


          Это сообщение было перенесено сюда или объединено из темы "Решение уравнений"
          user posted image
            Методы решения дифференциальных уравнений.

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

            ExpandedWrap disabled
                  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. Метод Эйлера.

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

            ExpandedWrap disabled
                  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]. Соотношения, описывающие
            данный метод, имеют вид:

            ExpandedWrap disabled
                  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 угла наклона касательной в точке:

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


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

            ExpandedWrap disabled
                  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. Метод Рунге-Кутта четвертого порядка.

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

            ExpandedWrap disabled
                  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. Метод Адамса-Бошфора.

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

            ExpandedWrap disabled
                  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):

            ExpandedWrap disabled
                  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, ...
            Итерационный процесс прекращается, когда
            ExpandedWrap disabled
                  | y^(i+1)_[m+1] - y^(i)_[m+1] | < Eps               (17)


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

            ExpandedWrap disabled
              {
                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 -
              Вычисление приближённого значения функций с помощью ряда Тейлора (Разложение функций в степенные ряды) (применяется также для нахождения пределов вообще)

              К примеру, 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).
              Для повторения вычисления этих функций смотрите прикреплённые вверху ссылки.

              Теперь надо собрать все эти функции в одном выражении, которое будет вычисляться в отдельной функции.
              ExpandedWrap disabled
                 { вычисление значения элемента с индексом [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.
              ExpandedWrap disabled
                function FunctionValue (x: single): single;
              На входе задаётся точка (число), в которой нужно вычислить значение функции.

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

              ExpandedWrap disabled
                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 -

              Прикреплённый файлПрикреплённый файлtailor.pas.zip (0.89 Кбайт, скачиваний: 455)
              1 пользователей читают эту тему (1 гостей и 0 скрытых пользователей)
              0 пользователей:


              Рейтинг@Mail.ru
              [ Script Execution time: 0,2298 ]   [ 15 queries used ]   [ Generated: 20.07.19, 11:42 GMT ]