Пусть дано дифференциальное уравнение

с начальным условием

.

Выбрав шаг h, положим . Приближенное решение заданного уравнения может быть полученно при помощи методов Рунге – Кутта, к которым относятся методы, обладающие следующими свойствами:

1.) Эти методы являются одношаговыми.

2.) Приближенные решения, найденные при помощи этих методов, совпадают с рядом Тейлора вплоть до членов порядка h^р, где степень р различна для различных методов и называется порядком метода.

3.) Эти методы не требуют вычисления производных от правой части дифференциального уравнения f(x,y), а требуют только вычисления самой функции.

Отметим, что для продвижения на один шаг (т.е. для вычисления у_n+1 по предыдущей точке (x_n,y_n) функцию f(x,y) придется вычислять в нескольких вспомогательных точках отрезка [x_n, x_n+1]. Однако это намного удобнее, чем вычислять значения производных.

Очевидно, что метод Эйлера, рассмотренный в параграфе 6.3, является одним из методов Рунге – Кутта первого порядка. Для продвижения на один шаг при помощи метода Эйлера требуется вычислить (см.(6.13)). Уточненный метод Эйлера, рассмотренный в пункте 6.4, является одним из методов Рунге – Кутта второго порядка. Для продвижения на один шаг при помощи уточненного метода Эйлера необходимо вычислить и (см. (6.14) – (6.16)).

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

(6.19)

Полагаем

. (6.20)

Выбор постоянных производим так, чтобы разложения (6.12) и (6.20) по степеням h совпадали до возможно более высоких степеней h при произвольной функции f(x,y) и произвольном шаге h (при разложении по степеням h используют формулу Тейлора для функции двух переменных).

Доказано, что при

(6.21)

разложения (6.12) и (6.20) по степеням h совпадают до членов порядка h⁴ включительно.

С учетом (6.21) формула (6.20) запишется в виде

, (6.22)

где вычисляются по формулам (6.19).

Таким образом, для продвижения на один шаг из точки необходимо по формулам (6.19) вычислить и по формуле (6.22) - .

Погрешность метода на одном шаге есть величина порядка h^s.

Рассмотренный метод, как уже указывалось, является методом Рунге – Кутта четвертого порядка. Этот метод применяется настолько широко, что в литературе он часто называется методом Рунге – Кутта без указаний на порядок.

Вычисления по методу Рунге – Кутта, необходимые для получения у₁ по известному у₀, занесены в табл. 6.1.

Таблица 6.1

n	x	y	k=hf(x,y)	Dyn
	x0	y0	k1(0)	Dy0=1/6(k1(0)+2 k2(0)+2 k3(0)+ k4(0))
	x0+h/2	y0+ k1(0)/2	k2(0)
	x0+h/2	y0+ k2(0)/2	k3(0)
	x0+h	y0+ k3(0)	k4(0)
	x1=x0+h	y1= y0+Dy0

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

. (6.23)

Величина не должна превышать нескольких сотых. В противном случае шаг h следует уменьшить.

Пример. Дано дифференциальное уравнение

с начальным условием у(1)=1.

При помощи метода Рунге – Кутта вычислить у(1,1). Шаг h принять равным 0,1. Вычисления занести в табл. 6.2.

Таблица 6.2

n	x	y	k=hf(x, y)
			0.2	0,209999
	1,05	1,1	0,20975
	1,05	1,104875	0,210237
	1,1	1,210237	0,220024
	1,1	1,209999

Итак, .

Функция у=х² является решением рассматриваемого уравнения. Таким образом, при х=1,1 точное значение решения следующее:

у(1,1)=1,21.

6.6. Экстраполяционный метод Адамса

Недостатком метода Рунге – Кутта является то, что для получения одного нового значения решения дифференциального уравнения правую часть уравнения приходится подсчитывать в нескольких точках. Если правая часть сложна, то требуется большой объем вычислительной работы. Рассмотрим многошаговые методы, применение которых на одном шаге не требует большого объема вычислительной работы. Пусть требуется найти решение дифференциального уравнения

с начальным условием

.

Предположим, что нам удалось каким-либо образом найти приближенное решение в точках . Обозначим

.

Большое количество различных формул численного решения заданного дифференциального уравнения можно получить следующим образом. Проинтегрируем обе части дифференциального уравнения от (i – любое целое число от 0 до m) до . Получим

или

(6.24)

Значения подынтегральной функции можно вычислить в (m+1)-й точке (). По вычисленным значениям строим интерполяционный многочлен степени m. Подставляя интерполяционный многочлен вместо f(x,y(x)) в (6.24), приближенно получим

(6.25)

Выбирая различные значения m, i и различные формы записи интерполяционного многочлена , получим различные формулы численного решения дифференциального уравнения.

Пусть i=0, т.е. нижний предел интегрирования в формуле (6.25) равен х_n. Равенство (6.25) запишется в виде

. (6.26)

Для записи интерполяционного многочлена воспользуемся второй интерполяционной формулой Ньютона

, (6.27)

где .

Подставляем (6.27) в (6.26):

, dx=hdq;

= q=0 при x=x_n; =

q=1 при x=x_n+1.

где

и т.д.

Таким образом

(6.28)

Эта формула называется экстраполяционной формулой Адамса (Дж. Адамс (1819 – 1892) – английский астроном, математик и механик) в связи с тем, что интерполяционный многочлен используется вне отрезка, содержащего узлы интерполяции.

В вычислительной практике чаще всего используется экстраполяционная формула Адамса при m=3, т.е. подынтегральная функция f(x,y(x)) в формуле (6.24) заменяется интерполяционным многочленом третьей степени. При построении интерполяционного многочлена используются четыре значения приближенного решения: y_n-3, y_n-2, y_n-1, y_n.

Экстраполяционная формула Адамса в этом случае имеет следующий вид:

(6.29)

Для начала вычислений по формуле (6.29) нужны четыре значения:

у₀, y₁, y₂, у₃. Значение у₀ задается начальным условием. Значения y₁, y₂, у₃ вычисляются каким-либо другим методом, например методом Рунге – Кутта.

Для применения формулы (6.29) необходимо составить таблицу конечных разностей функции f. При работе на ЭВМ оперирование таблицами конечных разностей вызывает затруднения. Преобразем формулу (6.29). конечные разности , , можно выразить через табличные значения функции f, а именно:

(6.30)

Подставляя (6.30) в (6.29), получаем экстраполяционную формулу Адамса в следующем виде:

. (6.31)

Оценим погрешность экстраполяционного метода Адамса (см. (6.29) и (6.31) на одном шаге. Для получения расчетных формул (6.29) и (6.31) использован интерполяционный многочлен третьей степени.

Следовательно, погрешность формул (6.29) и(6.31) можно записать в виде

.

Величина как погрешность второй интерполяционной формулы Ньютона пропорциональна h⁴. Следовательно, погрешность экстраполяционного метода Адамса на одном шаге есть величина порядка h⁵.

6.7. Метод Милна

Рассмотрим дифференциальное уравнение

с начальным условием

.

Выбрав шаг h, положим . Предположим, что нам удалось каким-либо образом найти приближенное решение в точках . Покажем, как по этим четырем значениям можно найти приближенное значение решения дифференциального уравнения в точке x_n+1.

Проинтегрируем обе части дифференциального уравнения от х_n-3 до х_n+1:

.

Отсюда

(6.32)

Для подынтегральной функции f(x,y(x)) построим интерполяционный многочлен третьей степени на отрезке по узлам . Воспользуемся первой интерполяционной формулой Ньютона

,

где

.

Подставляя интерполяционный многочлен вместо функции f в (6.32), приближенно получаем

= , dx=hdq;

q=0 при x=x_n-3; =

q=4 при x=x_n+1.

(6.33)

Преобразуем полученную формулу.

Выразим конечные разности через значения функции:

Подставляя полученные выражения для конечных разностей в формулу (6.33), получаем первую формулу Милна

. (6.34)

Для вывода второй формулы Милна интегрируем обе части исходного дифференциального уравнения от x_n-1 до x_n+1:

Отсюда

(6.35)

Для подынтеральной функции f(x,y(x)) построим интерполяционный многочлен третьей степени на отрезке . Воспользуемся первой интерполяционной формулой Ньютона

,

где

.

Подставляя интерполяционный многочлен вместо функции f в (6.35) приближенно получаем

= , dx=hdq;

q=0 при x=x_n-1; =

q=2 при x=x_n+1.

(6.36)

Преобразуем полученную формулу. Выразим конечные разности через значения функции:

.

Подставляя полученные выражения для конечных разностей в формулу (6.36), получим вторую формулу Милна

. (6.37)

Обозначим:

- значение у_к, найденные по первой формуле Милна;

.

Метод Милна применяется следующим образом.

1. Вычисляем первое приближение по первой формуле Милна:

.

2. По значению вычисляем

.

3. Находим второе приближение у_n+1 по второй формуле Милна:

.

Первая формула Милна служит «предсказывающей» формулой (формулой прогноза), вторая – «поправочной» формулой (формулой корректировки). Методы, в которых сначала находится предварительное значение функции по одной формуле, а затем это значение уточняется по другой формуле, объединяются под общим названием методов прогноза и коррекции.

Метод Милна, как видно из рассмотренного выше, - многошаговый метод. Для применения метода Милна необходимо найти первые четыре значения решения дифференциального уравнения: у₀, у₁, у₂, у₃, используя начальное условие и какой-либо метод, например, метод Рунге – Кутта.

Оценим погрешность метода Милна. Так же, как и для экстраполяционного метода Адамса, погрешность метода Милна на одном шаге есть величина порядка h^s.

На практике для оценки погрешности метода Милна ограничились разностями третьего порядка. Обозначим - погрешности первой и второй формул Милна соответственно. Учитывая отброшенные в интерполяционной формуле Ньютона разности четвертого порядка, с точностью до разностей пятого порядка будем иметь:

.

Отсюда, считая, что четвертые разности практически постоянны на отрезке длины 4h, получим

.

Так как

,

получаем

.

Отсюда

. (6.38)

6.8. Приближенное решение систем дифференциальных

уравнений и дифференциальных уравнений высших порядков

Рассмотрим систему трех дифференциальных уравнений:

(6.39)

где х – независимая переменная;

y, z, u – искомые функции.

Требуется найти решение системы, удовлетворяющее начальным условиям

. (6.40)

Обозначим

; . (6.41)

После введения новых обозначений система (6.39) запишется в виде

, (6.42)

а начальные условия (6.40) – в виде

. (6.43)

Все рассмотренные выше численные методы решения задачи Коши для дифференциального уравнения первого порядка переносятся без изменений на случай системы уравнений. Формальное отличие состоит лишь в том, что в соответствующих соотношениях вместо скалярных величин у, у₀, f(x,y) участвуют векторные величины Y, Y₀, F(x,Y). Например, расчетная формула метода Эйлера для решения заданной системы (6.42) с начальным условием (6.43) запишется в виде

. (6.44)

где .

В развернутом виде формула (6.44) запишется следующим образом:

(6.45)

аналогично поступают и при решении системы другими методами. Например, расчетная формула метода Рунге – Кутта имеет следующий вид:

) ,

где

В развернутом виде каждое из этих соотношений запишется в виде трех соотношений для скалярных величин.

Таким образом между решением одного дифференциального уравнения первого порядка и решением системы дифференциальных уравнений имеется лишь часто количественное различие – число требуемых вычислений растет пропорционально числу неизвестных функций.

Пример. Применяя метод Эйлера, найти приближенное решение системы дифференциальных уравнений

с начальными условиями

у(0)=0; z(0)=-1

на отрезке [0, 0.3] с шагом y=0.1. Вычисления оформить в виде табл. 6.3.

Таблица 6.3

		-1		0.1
	0.1	0.1 -1	1.2 0.4	0.12 0.04
	0.2	0.22 -0.96	1.4 0.8	0.14 0.08
	0.3	0.36 -0.88

Приближенное решение системы в последовательных точках находим по формуле (6.44):

;

, т.е.

, т.е.

, т.е.

Отметим, что решением данной системы являются функции . Сравним полученные приближения Y₁, Y₂, Y₃ с точными:

; .

Рассмотрим дифференциальное уравнение n-го порядка

с начальными условиями

; ; .

Обозначим

; ; …; .

В этом случае вместо заданного дифференциального уравнения можно рассматривать эквивалентную систему дифференциальных уравнений:

……………….

с начальными условиями

; ; ; .

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

6.9. Краевые задачи для обыкновенных дифференциальных

уравнений. Постановка задачи

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

Пусть дано дифференциальное уравнение второго порядка

. (6.46)

Двухточечная краевая задача для уравнения (6.46) ставится следующим образом. Найти функцию y=f(x), которая внутри отрезка [a,b] удовлетворяет уравнению (6.46), а на концах отрезка – краевые условия

;

.

Приведем некоторые примеры краевых задач.

1. Найти функцию y=y(x), удовлетворяющую внутри отрезка [a,b] уравнению

и принимающую на концах отрезка заданные значения

у(а)=А, y(b)=B.

Геометрически это означает, что требуется найти интегральную кривую, проходящую через данные точки (а, А) и (b, B).

2. Найти функцию y=y(x), удовлетворяющую внутри отрезка [a,b] уравнению

,

причем производная на концах отрезка должна принимать заданные значения

, .

Геометрически это означает, что требуется найти интегральную кривую, пересекающую прямые х=а и х=b под заданными углами и соответственно, такими, при которых , .

Если дифференциальное уравнение и краевые условия линейны, краевая задача называется линейной краевой задачей. В этом случае дифференциальное уравнение и краевые условия записываются так:

где Р(х), q(x), f(x) – известные непрерывные на отрезке [a,b] функции;

- заданные постоянные, причем и .

6.10. Решение краевых задач

методом конечных разностей

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

(6.47)

с условиями

. (6.48)

Разобъем отрезок [a,b] на n равных частей точками

.

Точки x_i будем называть узлами.

Рассмотрим дифференциальное уравнение (6.47) во внутренних узлах отрезка, т.е. в точках . Полагаем в (6.47) х= x_i. Тогда

. (6.49)

Обозначим получаемые в результате приближенные значения искомой функции y(x) в узлах x_i через y_i. Введем также обозначения:

.

В каждом внутреннем узле производную можно приближенно заменить, используя одну из формул:

;

(6.50)

.

Аналогично производную можно приближенно заменить, используя одну из формул:

(6.51)

.

Для получения формул (6.50) и (6.51) можно воспользоваться интерполированием и численным дифференцированием.

Для точки х₀=а производную приближенно заменяем, используя первую формулу из (6.50):

, (6.52)

а для точки х_n=b производную заменяем, используя вторую формулу из (6.50):

. (6.53)

Отметим, что заменяем именно по этим формулам, так как это крайние точки отрезка. Для внутренних узлов отрезка [a,b] можно использовать любую формулу из (6.50) и (6.51). На практике чаще используется третьи формулы из (6.50) и (6.51), как более точные.

Подставляя и из этих формул в уравнения (6.49), получим (n-1) уравнение. Подставляя (6.52) и (6.53) в краевые условия (6.48), получим еще два уравнения.

Итак, после подстановок приходим к линейной системе из (n++1)-го уравнения с (n+1)-й неизвестной:

(6.54)

.

Решив систему, если это возможно, получим таблицу приближенных значений искомой функции. При большом n непосредственное решение системы (6.54) становится затруднительным. Для решения систем уравнений вида (6.54) разработан специальный метод, получивший название метода прогонки.

6.11. Метод прогонки

Рассмотрим линейную систему уравнений (6.54). Преобразуя первые (n-1) уравнения:

или

.

Отсюда

, (6.55)

где

(6.56)

Предположим, что с помощью полной системы (6.55) из i-го уравнения системы исключен член, содержащий . Тогда это уравнение можно записать в виде

(6.57)

Аналогично

.

Подставляя выражение для в i-е уравнение системы (6.55), получим

.

Отсюда

. (6.58)

Сравнивая (6.57) с (6.58), получим

. (6.59)

Получены рекуррентные соотношения для определения и . Для начала подсчетов необходимо знать и . Рассмотрим систему уравнений:

(получено из (6.55) при i=1);

(получено из первого краевого условия).

Находим у₀ из второго уравнения системы:

(6.60)

и, подставляя его в первое уравнение системы, получаем

.

Сравнивая полученное выражение с выражением

,

получаем

; (6.61)

Итак, зная и , по формулам (6.59) можно последовательно вычислить коэффициенты и для . Вычисление коэффициентов и называется прямым ходом.

Определим y_n из системы уравнений:

(получено из (6.57) при i=n-1);

(получено из второго краевого условия);

. (6.62)

Теперь, используя формулы (6.57) и (6.60), последовательно находим:

, , …, ;

……………………………

.

Вычисление y_i называется обратным ходом.

Вычисления удобно оформить в виде табл. 6.4.

Таблица 6.4

					Прямой ход	Обратный ход
…	…	…	…	…	…	…	…
n-1
n

Глава 7. ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ В ЧАСТНЫХ ПРОИЗВОДНЫХ

7.1. Классификация дифференциальных уравнений

в частных производных

Дифференциальные уравнения в частных производных имеют широкое применение в различных областях науки и техники. Получить их решение в явном виде удается только в самых простейших случаях. В связи с этим важное значение приобретают приближенные методы решения различных задач для дифференциальных уравнений в частных производных. На практике чаще всего используются линейные дифференциальные уравнения второго порядка.

Пусть D – некоторая область изменения независимых переменных х и у, ограниченная контуром Г. в области задано линейное дифференциальное уравнение второго порядка для функции U(x, y), если для любой точки из области D имеет место соотношение

. (7.1)

где - коэффициенты;

f(x,y) – свободный член уравнения.

Если коэффициенты a, b, c, d, e, g не зависят от х и у, то уравнение (7.1) является линейным дифференциальным уравнением с постоянными коэффициентами.

Обозначим

.

В зависимости от знака функции линейное дифференциальное уравнение (7.1) относится в данной области D к одному из следующих типов:

К эллиптическому - <0;

К параболическому - =0;

К гиперболическому - >0.

В дальнейших параграфах будут рассмотрены примеры каждого из трех типов уравнений и методы их приближенного решения.

7.2. Метод сеток решения краевых задач

для дифференциальных уравнений эллиптического типа

Рассмотрим дифференциальное уравнение эллиптического типа

(7.2)

Это уравнение называют уравнением Пуассона (С. Пуассон (1781 – 1840) – французкий механик, физик и математик).

Первая краевая задача, или задача Дирихле для уравнения Пуассона, ставится следующим образом. Найти функцию u=u(x,y), удовлетворяющую внутри некоторой области D уравнению (7.2), а на границе Г – условию

, (7.3)

где - заданная непрерывная функция.

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

Проводим на плоскости х0у два семейства параллельных прямых

; .

Точки пересечения этих прямых назовем узлами. Узлы имеют координаты:

, .

Условимся на рисунках данной главы указывать только индексы узла. Например, узел (x_i, y_k) обозначается через (i, k).

Два узла называются соседними, если они удалены друг от друга в направлении оси 0х или 0у на расстояние, равное шагу сетки h или l соответственно. Таким образом, расположение узла (x_i, y_k) и его соседних узлов выглядит следующим образом (7.1).

(i, k+1)

	l

h

(i-1,k) (i,k) (i+1,k)

(i,k-1)

Рис. 7.1

Будем рассматривать только те узлы, которые принадлежат области . Разобъем узлы области на внутренние и граничные. Некоторый узел (x_i, y_k) будем считать внутренним, если он сам или четыре соседних узла принадлежат области . Узлы области , у которых хотя бы один соседний узел не принадлежит области , называется граничным. Множество внутренних узлов обозначим через , множество граничных узлов – через Г*.

На рис. 7.2 внутренние узлы помечены знаком «о», граничные – знаком «*»

у

0 х

Рис. 7.2

Рассмотрим дифференциальное уравнение (7.2) во всех внутренних узлах:

. (7.4)

Заменим частные производные второго порядка разностными отношениями:

(7.5)

где через U_ik обозначено приближенное сеточное значение решения уравнения (7.2) в узле (x_i, y_k).

Подставляя (7.5) в (7.4), получим для каждого внутреннего узла приближенное равенство

. (7.6)

Рассмотрим граничные узлы, т.е. узлы (x_i, y_k) Г*. Известно, что по условию (7.3). Поэтому, если узел попадает на границу Г области D, значение U_ik можем вычислить, используя это условие. Однако, как правило, граничные узлы не будут попадать на границу Г. Пусть узел (x_i, y_k) не попадает на границу Г. В этом случае значение U_ik положим равным значению функции в точке границы Г, ближайшей к этому узлу.

Иногда значение U_ik полагают равным значению функции в точке Г, ближайшей к рассматриваемому гра ничному узлу в направлении оси 0х (рис. 7.3):

U_ik=U(B) (x_i, y_k) Г*. (7.7)

Объединяя равенства (7.6) и (7.7), получим неоднородную линейную систему, в которой число уравнений равно числу неизвестных. Доказанно, что эта система имеет решение и притом единственное. Решив ее, получим приближенные значения искомого решения в узлах сетки.

Рассмотрим частный случай уравнения (7.2), который получается при f(x,y)=0

. (7.8)

Это уравнение называют уравнением Лапласа (П. Лаплас (1749 – 1827) – французский математик, физик и астроном).

Конечно – разностные уравнения (7.6) для уравнения Лапласа запишутся в виде

. (7.9)

Наиболее простой вид эти уравнения имеют при h=l:

или . (7.10)

у

М В

0 ч

Рис. 7.3

Объединяя (7.7) и (7.10), получим линейную систему уравнений. Решив систему, получаем приближенные значения решения уравнения (7.8), удовлетворяющего условию (7.3), в узлах сеточной области.

7.3. Погрешность аппроксимации дифференциальных уравнений

эллиптпческого типа конечно – разностыми уравнениями

При применении метода сеток для решения краевых задач дифференциальные уравнения заменяются конечно – разностными уравнениями. Выясним, чему равна погрешность такой замены. Рассмотрим, например, погрешность, получаемую в результате замены дифференциального уравнения (7.8) конечно – разностным уравнением (7.10).

Используя формулу Тейлора, значения функции , , , запишем в виде

(7.11)

где

, , ,

.

Складывая равенства (7.11), получим

Отсюда

,

где

.

Если ввести обозначение

,

то

.

Таким образом, погрешность замены дифференциального уравнения (7.8) конечно – разностным уравнением (7.10) имеет порядок h². Такой же порядок имеет погрешность замены дифференциального уравнения (7.8) конечно – разностным уравнением (7.9), а также погрешность замены дифференциального уравнения (7.2) конечно – разностным уравнением (7.6).

7.4. Аппроксимация граничных условий

В параграфе 7.2. указывалось, что в граничных узлах значения функции можно заменить значением функции в точке границы Г, ближайшей к рассматриваемому граничному узлу в направлении оси 0х.

Очевидно, что погрешность этой замены имеет порядок h, так как для вычисления U(B) был использован интерполяционный многочлен нулевой степени. Точность вычисления U(B) можно повысить, если для вычисления U(B) использовать интерполяционный многчлен первой степени. Для построения такого многочлена, кроме точки М, используем внутренний узел А, ближайший к узлу В по направлению оси 0х (рис. 7.4):