Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Методы решения дифференциальных уравнений
Данная статья посвящена численным методам решения дифференциальных уравнений. Мы будем рассматривать методы решения одного обыкновенного дифференциального уравнения первого порядка с одним начальным условием:
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.
Метод Эйлера
Этот метод, один из самых старых и широко известных, описывается формулой:
y[m+1] = y[m] + h*f(x[m],y[m]). (3)
Найденное по формуле (3) решение согласуется с разложением в ряд Тейлора вплоть до членов порядка h, т.е. данный метод является методом Рунге-Кутта первого порядка.
Этот метод имеет довольно большую ошибку приближения, кроме того, он очень часто оказывается неустойчивым - изначально малая ошибка (происходящая от приближения, округления или исходных данных) увеличивается с ростом x.
Для вычисления y[m+1] метод Эйлера использует наклон касательной только в точке x[m],y[m]. Этот метод можно усовершенствовать множеством различных способов.
Рассмотрим два из них.
Исправленный метод Эйлера
В исправленном методе Эйлера мы находим средний tg угла наклона касательной для двух точек: x[m],y[m] и x[m+1],y[m]+h*y’[m]. Соотношения, описывающие данный метод, имеют вид:
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, являясь, т.о. методом Рунге-Кутта второго порядка.
Модифицированный метод Эйлера
В данном методе мы находим 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, и также является методом Рунге-Кутта второго порядка.
Метод Рунге-Кутта четвертого порядка
Данный метод является одним из самых употребительных методов интегрирования дифференциальных уравнений. Этот метод применяется настолько широко, что в литературе его просто называют “методом Рунге-Кутта” без всяких указаний на его тип или порядок. Этот классический метод описывается системой следующих пяти соотношений:
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], а затем используется тот или иной метод его “корректировки”.
Среди множества возможных формул прогноза и коррекции выберем по одному примеру, применимому ко многим практическим задачам.
Метод Адамса-Бошфора
Для прогноза используем формулу второго порядка
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.
{$E+,N+,X+}
const
eps = 1e-5;
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.