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