Численные методы , Уравнения, СЛАУ, пределы, ряды

[Romtek]

Решение уравнений

{Автор: 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]

Численное решение уравнений итерационными методами.

Автор:

  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 - 06.11.06, 12:09

[e-moe]

Численное решение систем линейных алгебраических уравнений

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

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. Сравнение методов.

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

• Время вычислений пропорционально n² на итерацию, в то время как для метода исключений время пропорционально n³. Если для решения системы требуется меньше чем 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.

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

[Romtek]

Методы решения дифференциальных уравнений.

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

    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 - 30.06.05, 09:09

[Romtek]

Вычисление приближённого значения функций с помощью ряда Тейлора (Разложение функций в степенные ряды) (применяется также для нахождения пределов вообще)

К примеру, sin = x -(x³/3!)+(x⁵/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)ⁿ * x²ⁿ⁺¹ / (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 - 30.06.05, 09:14

Прикреплённый файлtailor.pas.zip (0.89 Кбайт, скачиваний: 618)