Тема 5. Интегрирование дифференциальных уравнений.

1. Постановка задачи.

Дано уравнение F(x, y, y', y''…y⁽ⁿ⁾)=0 (1). Здесь x – аргумент, y – функция. Общий интеграл (общее решение) записывается в виде Φ(x, y, c₁, c₂…c_n)=0 (2). Здесь c_i – произвольные постоянные. Решение (2) при фиксированных c_i называется частным решением. Для выделения из общего решения частного необходимо задать n условий, называемых начальными. Пусть начальные условия заданы в виде y^(k)(x₀)=y_0k, где k=0, 1…n-1 (3). Задачей Коши называется задача интегрирования (1) с начальными условиями (3). Иначе задача Коши называется задачей с начальными условиями. Численное решение задачи Коши – это отыскание для x₀, x₁…x_n значений f(x₀), f(x₁)…f(x_n), т. е. построение таблицы y(x). Численные методы интегрирования дифференциальных уравнений делятся на шаговые (Эйлера, Рунге-Кутта, Гилла) и разностные (Адамса, Милна, Штернера, Гаусса) методы. Обычно уравнение (1) разрешают относительно старшей производной, т. е. y⁽ⁿ⁾=f(x, y, y',y''…y^(n-1)) (4). Его общий интеграл запишется в виде y=φ(x, c₁, c₂…c_n) (5). Таким образом, займемся решением задачи Коши для уравнения (4) с начальными условиями (3).

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

Пусть дано дифференциальное уравнение первого порядка. =y'=f(x, y(x)) (6) и начальное условие y(x₀)=y₀ (7). Решение уравнения (5) в конечном итоге сводится к интегрированию обыкновенного дифференциального уравнения вида (6) с начальным условием (7). Существование и единственность решения дается специальной теоремой из математического анализа. Будем предполагать, то ее условия выполняются, т. е. решение существует. Построим систему равноотстоящих узлов x₀, x₁…x_n, где x_i=x₀+ih, h=const>0. Проинтегрируем (6) на отрезке [x_i, x_i+1].

'dx= (x, y(x))dxy(x_i+1)=y(x_i)+ (x, y(x))dx.

Если в интеграле принять f(x, y(x))=const=f(x_i, y(x_i)), то по квадратурной формуле левых прямоугольников получим y(x_i+1)≈y(x_i)+(x_i+1-x_i)•f(x_i, y(x_i)); y(x_i+1)≈y(x_i)+h•f(x_i, y(x_i)).

Это шаговая формула, т. е. для вычисления последующего значения необходимо знать предыдущее. Перейдем к сеточному представлению функции (т. е. приближенному решению) заменив y_i≈y(x_i). Получим y_i+1=y_i+Δy_i; Δy_i=h•f(x_i, y_i) x_i=x₀+ih; i=0, 1…N (8).

Это формула Эйлера. Геометрическая интерпретация метода довольно проста. В точке M(x₀, y₀) к отыскиваемой функции строится касательная. Затем ищется точка M₁, как пересечение касательной и прямой x=x₁. Пусть O₁ – точка пересечения прямых x=x₁ и y=y₀. Рассмотрим треугольник MO₁M₁. O₁M₁=MO₁•tg(O₁MM₁)=h•y'(x₀)=h•f(x₀, y₀); - = -y₀=h•f(x₀, y₀); =y₀+h•f(x₀, y₀)=y₁. Следующий шаг – построение касательной в точке M₁ к кривой из семейства решений и отыскание точки пересечения этой касательной с прямой x=x₂. Затем – построение касательной в точке M₂ и так далее. Таким образом, в зависимости от начальных условий, получим частное решение семейства. MM₁M₂… называется ломаной Эйлера. Пусть M₁' – точка на истинной кривой решения. Тогда M₁M₁' – погрешность метода (совпадает с шагом метода).

3. Сходимость и устойчивость метода Эйлера.

Пусть есть уравнение (6) y'=f(x, y) и начальное условие y(x₀)=y₀. Если x_i=x₀+ih, то y_i+1=y_i+h•f(x_i, y_i) (8). Пусть y(x_i) – точное решение (6) в точке x_i, а y_i – точное решение уравнения (8). Погрешность метода на i-том шаге ε_i=y_i-y(x_i). Метод называют сходящимся, если ε_i→0 при h→0, т. е. сходимость – уменьшение погрешности с уменьшением шага. Доказана теорема, что h→0 на любом конечном отрезке метод Эйлера сходится, т. к. ε_i< •(exp(k(x_i-x₀))-1). Здесь M и k – некоторые оценочные константы. Неустойчивостью называют накопление погрешности вычислений в результате неточности вычислений на каждом шаге. Метод Эйлера устойчив. Доказано что если на первом шаге погрешность вычислений = -y₀, то на i-том шаге | |≤exp(k(x_i-x₀))•| |, т. е. погрешность имеет тот же порядок малости, что и шаг (ε_i≈h). Таким образом, Метод Эйлера имеет первый класс точности.

4. Усовершенствования метода Эйлера.

Легко видеть, что при всей его сходимости и устойчивости обладает слишком большой погрешностью. В частности, если производная функции знакопостоянна и положительна, то решение по методу Эйлера будет все время ниже истинного решения. Очевидно, можно уменьшать шаг h до тех пор, пока погрешность не станет меньше заданной. Но это приведет к столь малому шагу, что очень сильно возрастет объем вычислений, и метод потеряет свое основное преимущество – простоту. Этот недостаток метода связан с тем, что мы вычисляем значение функции только в текущей точке и не пытаемся предугадать, характер изменения функции в будущей точке.

4.1. Алгоритм Эйлера-Коши.

Пусть есть уравнение (6) y'=f(x, y) и начальное условие y(x₀)=y₀ (*) Построим численное решение на системе узлов x₀, x₁…x_n ,где x_i+1=x₀+ih. Проинтегрируем (*) на отрезке [x_i, x_i+1] 'dx= (x, y(x))dx. По формуле трапеций, y(x_i+h)-y(x_i)≈ •(f(x_i, y(x_i))+f(x_i+1, y(x_i+1))); y(x_i+h)≈y(x_i)+ •(f(x_i, y(x_i))+f(x_i+1, y(x_i+1))) (**)

Уравнение (**) можно рассматривать как уравнение, описывающее y(x_i+1) неявно, т. е. y(x_i+1)≈B•y(x_i+1).

Поскольку выразить y(x_i+1) удается крайне редко, то получим его по методу Эйлера и подставим в (**) y(x_i+h)≈y(x_i)+h•f(x_i, y(x_i)); y(x_i+h)≈y(x_i)+ •(f(x_i, y(x_i))+f(x_i+1, y(x_i)+h•f(x_i, y(x_i)))).

Перейдем к сеточному представлению, обозначив x_i+1=x_i+h; y(x_i)=y_i; y(x_i+1)=y_i+1.

Таким образом, формулы Эйлера-Коши имеют вид

y_i+1*=y_i+h•f(x_i, y_i); y_i+1=y_i+ •(f(x_i, y_i)+f(x_i+1, y_i+1*)).

Этот метод имеет второй класс точности, т. е. ε_i~h². Геометрический смысл метода состоит в следующем. Пусть A – точка с координатами (x₀, y₀) (начальная точка). В этой точке к кривой из семейства решений проводится касательная. Точка B – точка пересечения этой касательной с прямой x=x₁. В точке B проводится новая касательная к кривой из семейства решений. Если α – угол с положительным направлением оси Ox первой касательной, а β – угол с положительным направлением оси Ox второй касательной, то следующая точка будет пересечением прямой x=x₁ и прямой из точки B, угол с положительным направлением оси Ox которой γ удовлетворяет условию tg(γ)= . Пусть С – точка пересечения этой прямой с прямой x=x₁, а точка D – точка пересечения прямой y=y₀ и прямой x=x₁. Из треугольника ACD следует, что CD=AD•tg(γ)=h• = •(y'(x₀)+y'(x₁))= •(f(x₀, y₀)+f(x₁, y₁*)). Т. к. CD=y_C-y₀, то y_C=y₀+ •(f(x₀, y₀)+f(x₁, y₁*))=y₁.

4.2. Алгоритм модифицированного метода Эйлера.

Пусть есть уравнение (6) y'=f(x, y) и начальное условие y(x₀)=y₀ (*) Построим численное решение на системе узлов x₀, x₁…x_n ,где x_i+1=x₀+ih. Проинтегрируем (*) на отрезке [x_i, x_i+1] 'dx= (x, y(x))dx. По формуле центральных прямоугольников, y(x_i+h)-y(x_i)= (x, y(x))dx≈h•f(x_i+ , y(x_i+ )); y(x_i+h)≈y(x_i)+h•f(x_i+ , y(x_i+ )). После интегрирования появится неизвестная y(x_i+ ). Ее найдем методом Эйлера с шагом . y(x_i+ )≈y(x_i)+ •f(x_i, y(x_i)); y(x_i+h)≈y(x_i)+h•f(x_i+ , y(x_i)+ •f(x_i, y(x_i))). Перейдем к сеточному представлению, обозначив x_i+h/2=x_i+1/2; y(x_i+ )=y_i+1/2; y(x_i+h)=y_i+1. Таким образом, формулы модифицированного метода Эйлера имеют вид y_i+1/2=y_i+ •f(x_i, y_i); y_i+1=y_i+h•f(x_i+1/2, y_i+1/2). Этот метод также имеет второй класс точности. Геометрический смысл метода состоит в следующем. Пусть M – точка с координатами (x_i, y_i) (текущая точка). В этой точке к кривой из семейства решений проводится касательная. Пусть эта касательная имеет угол α_i с положительным направлением оси Ox. Из точки пересечения этой касательной и прямой x=x_i+1/2 проводится новая касательная, которая имеет угол α_i+1/2 с положительным направлением оси Ox. Новая точка M₁ будет точкой пересечения прямой x=x₁ с прямой, имеющей угол α_i+1/2 с положительным направлением оси Ox, проведенной из точки M,

4.3. Усовершенствованный метод Эйлера-Коши с последующей итерационной обработкой.

Метод более точен, чем метод Эйлера-Коши. На каждом шаге производится итерационная обработка каждого найденного значения y_i. В начале выбирается грубое приближение y_i⁽⁰⁾=y_i+h•f(x_i, y_i). Далее строится итерационный процесс y_i+1^(k)=y_i+ •(f(x_i, y_i)+f(x_i+1, y_i+1^(k-1))). Итерации продолжаются до тех пор, пока два приближения y_i+1^(k) и y_i+1^(k+1) не совпадут в интересующем вычислителя количестве знаков. После этого полагают y_i+1=y_i+1^(k).

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

Пусть имеется система y'=f(x, y, z),z'=g(x, y, z) с начальным условием y(x₀)=y₀, z(x₀)=z₀.

Перейдем к сеточному представлению, обозначив y(x_i)≈y_i; z(x_i)≈z_i. Тогда y_i+1=y_i+Δy_i; z_i+1=z_i+Δz_i; Δy_i=h•f₁(x_i, y_i, z_i); Δz_i=h•f₂(x_i, y_i, z_i).

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

Одним из методов повышенной точности является метод Рунге-Кутта. Этот метод основан на прогнозе и коррекции тангенса угла наклона касательной.

Пусть есть уравнение (6) y'=f(x, y), начальное условие y(x₀)=y₀ и отрезок интегрирования, который можно разбить на n частей, получив шаг h= . Общая идея шаговых методов состоит в том, что значение в следующей точке находится исходя из значения в предыдущей точке, т. е. y_i+1, как и в методе Эйлера, определяется по формуле y_i+1=y_i+Δy_i. Разложим y(x_i+h) в ряд Тейлора в окрестности точки x_i. y(x_i+h)=y(x_i)+h•y'(x_i)+ •y''(x_i)+ •y'''(x_i)+ •y⁽⁴⁾(x_i)+…

Ограничимся членами до h⁴ и вычислим Δy_i.

Δy_i=y(x_i+h)-y(x_i)=h•y'(x_i)+ •y''(x_i)+ •y'''(x_i)+ •y⁽⁴⁾(x_i), где y''(x_i), y'''(x_i), y⁽⁴⁾(x_i) находят дифференцированием уравнения (1).

Идея Рунге состоит в том, чтобы искать решение Δy_i в виде Δ_i=a₁k₁⁽ⁱ⁾+a₂k₂⁽ⁱ⁾+a₃k₃⁽ⁱ⁾+a₄k₄⁽ⁱ⁾+…, где k₁⁽ⁱ⁾=h•f(x_i, y_i); k₂⁽ⁱ⁾=h•f(x_i+b₁•h, y_i+c₁k₁⁽ⁱ⁾); k₃⁽ⁱ⁾=h•f(x_i+b₂•h, y_i+c₂k₂⁽ⁱ⁾); k₄⁽ⁱ⁾=h•f(x_i+b₃•h, y_i+c₃k₃⁽ⁱ⁾); …

При этом Δ_i должно совпадать с точностью до членов малости заданной степени h с остатком ряда Тейлора Δy_i=h•y'(x_i)+ •y''(x_i)+…

Рассмотрим идею метода на примере метода второго порядка, т. е. Δ_i=h•y'(x_i)+ •y''(x_i), а членами более высоких порядков малости пренебрежем.

Δ_iТейл=h•y'(x_i)+ •y''(x_i)+0(h³); Δ_i=a₁k₁⁽ⁱ⁾+a₂k₂⁽ⁱ⁾ (*); k₁⁽ⁱ⁾=h•f(x_i, y_i); k₂⁽ⁱ⁾=h•f(x_i+b₁•h, y_i+c₁k₁⁽ⁱ⁾). Индекс (i) у коэффициентов можно опустить.

Разложим k₂ в ряд Тейлора с точностью до второго порядка малости.

k₂=h•f_i+h²•b₁ +h•c₁k₁ +0(h³); k₁=h•f_i; k₂=h•f_i+h²•b₁ +h²•c₁f_i .

Подставим k₁ и k₂ в (*).

Δ_i=h•a₁f_i+h•a₂f_i+h²•a₂b₁ +h²•a₂c₁f_i ; Δ_iТейл=h•y_i'+ •y_i''=h•f_i+ • + •f_i .

Сравнивая Δ_i и Δ_iТейл, получим a₁+a₂=1; a₂b₁= ; a₂c₁= . Эта система имеет бесконечное множество решений. Если положить a₁=a₂= , то b₁=c₁=1. Тогда k₁=h•f(x_i, y_i) – формула Эйлера. k₂=h•f(x_i+h; y_i+k₁); Δ= •(f(x_i, y_i)+f(x_i+h; y_i+k₁)); y_i+1=y_i+Δ – формула Эйлера-Коши. Можно взять и другие значения коэффициентов a₁ и a₂, получив другие значения k₁ и k₂, но можно строго показать, что такие значения коэффициентов являются оптимальными.

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

Если Δ_i=a₁k₁+a₂k₂+a₃k₃+a₄k₄, то k₁=h•f(x_i, y_i); k₂=h•f(x_i+ , y_i+ ); k₃=h•f(x_i+ , y_i+ ); k₄=h•f(x_i+h, y_i+k₃), т. е. b₁=b₂=c₁=c₂= ; b₃=c₃=1.

Значения коэффициентов a₁=a₄= ; a₂=a₃= .

Таким образом, рабочие формулы метода k₁⁽ⁱ⁾=h•f(x_i, y_i); k₂⁽ⁱ⁾=h•f(x_i+ , y_i+ ); k₃⁽ⁱ⁾=h•f(x_i+ , y_i+ ); k₄⁽ⁱ⁾=h•f(x_i+h, y_i+k₃⁽ⁱ⁾); y_i+1=y_i+ (k₁⁽ⁱ⁾+2k₂⁽ⁱ⁾+2k₃⁽ⁱ⁾+k₄⁽ⁱ⁾). Метод имеет четвертый класс точности. Недостаток метода в том, что для вычисления каждого значения y_i+1 требуется четыре вычисления правой части дифференциального уравнения, поэтому если у уравнения правая часть сложна, то применяют более эффективные методы.

Метод можно применять и для систем дифференциальных уравнений. Пусть имеется система y'=f(x, y, z), z'=g(x, y, z) с начальным условием y(x₀)=y₀, z(x₀)=z₀. Определим значения Δy_i, Δz_i. Δy_i= (k₁⁽ⁱ⁾+2k₂⁽ⁱ⁾+2k₃⁽ⁱ⁾+k₄⁽ⁱ⁾); Δz_i= (e₁⁽ⁱ⁾+2e₂⁽ⁱ⁾+2e₃⁽ⁱ⁾+e₄⁽ⁱ⁾); k₁⁽ⁱ⁾=h•f(x_i, y_i); e₁⁽ⁱ⁾=h•g(x_i, y_i); k₂⁽ⁱ⁾=h•f(x_i+ , y_i+ ); e₂⁽ⁱ⁾=h•g(x_i+ , y_i+ ); k₃⁽ⁱ⁾=h•f(x_i+ , y_i+ ); e₃⁽ⁱ⁾=h•g(x_i+ , y_i+ ); k₄⁽ⁱ⁾=h•f(x_i+h, y_i+k₃⁽ⁱ⁾); e₄⁽ⁱ⁾=h•g(x_i+h, y_i+e₃⁽ⁱ⁾). y_i+1=y_i+Δy_i; z_i+1=z_i+Δz_i.

6. Экстраполяционный метод Адамса (Разностный метод Адамса).

При решении дифференциальных уравнений методом Рунге-Кутта требуется многократные вычисления правой части для каждого y_i+1. Если правая часть дифференциального уравнения сложна, то применение метода Рунге-Кутта неоправданно. Метод Адамса не требует многократного вычисления правой части дифференциального уравнения. Пусть имеется уравнение y'=f(x, y(x)) с начальным условием y(x₀)=y₀, которое требуется проинтегрировать на отрезке [a, b]. Построим систему равноотстоящих узлов, разбив [a, b] на n частей. Тогда x_i=x₀+ih. Проинтегрируем заданное дифференциальное уравнение на промежутке [x_i, x_i+1]. y_i+1=y_i+ 'dx или Δy_i= 'dx. Заменим y' интерполяционным полиномом Ньютона в конце таблицы, ограничившись тремя конечными разностями.

y'(t)=y'_i+tΔy'_i-1+ •Δ²y'_i-2+ •Δ³y'_i-3, где t= .

Если x(x_i, x_i+1), то t(0, 1); dx=hdt. Подставим y'(t) в выражение для Δy_i, получим Δy_i=h• y'_i+tΔy'_i-1+ •Δ²y'_i-2+ •Δ³y'_i-3)dx=h•(y'_i+ •Δy'_i-1+ •Δ²y'_i-2+ •Δ³y'_i-3) – формула Адамса. Эту формулу можно преобразовать.

Δy_i= •(24y'_i+12•Δy'_i-1+10•Δ²y'_i-2+9•Δ³y'_i-3).

Т. к. Δy'_i-1=y'_i-y'_i-1; Δ²y'_i-2=y'_i-2y'_i-1+y'_i-2; Δ³y'_i-3=y'_i-3y'_i-1+3y'_i-2-y'_i-3, рабочая формула метода Адамса имеет вид y_i+1=y_i+ •(55y'_i-59y'_i-1+37y'_i-2-9y'_i-3).

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