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

 

Уравнение

F(x, у, у', у", ... , у(n))= 0,

связывающее неизвестную функцию у(х),независимую переменную х и производные у'(x), у"(x),..., у(n)(x) неиз­вестной функции, называется обыкновенным дифферен­циальным уравнением. Порядок п старшей производной называется порядком дифференциального уравнения.

В задаче Коши для дифференциального уравнения п-го порядкаискомая функция у(х),кроме самого дифферен­циального уравнения, должна удовлетворять начальным условиям

у(х0) = α0, у'(х0)= α1, у"(х0) = α2,..., у(n-1)(x0) = αn-1.

Методы решения задач для дифференциальных урав­нений можно разбить на три типа: точные, приближен­ные и численные.

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

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

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

7.1. Задача Коши

Рассмотрим задачу Коши для обыкновенного диффе­ренциального уравнения первого порядка

y'(х) = f(x, у), (7.1)

y(x0) = y0 (7.2)

Требуется найти функцию у = у(х), которая удовлет­воряет уравнению (7.1) на интервале (х0, X) и начально­му условию (7.2) в точке х0.

Приведем без доказательства теорему существования и единственности решения задачи Коши.

Теорема 7.1. Пусть в области R{(x, у), |х - х0| £ а, |у - у0| £ b} функция f(x, у) непрерывна. Тогда на некото­ром отрезке |х - x0| £ d существует решение уравнения (7.1), удовлетворяющее условию (7.2).

Если в области R функция f(x, у) удовлетворяет усло­вию Липшица

| f(x, у1) - f(x, y2)| £ k | y1 - у2|, k > 0, k = const,

то указанное решение единственно.

Произведем разбиение отрезка [х0, X] на п частей:

xi = x0 + ih, . (7.3)

Найдем приближенные значения решения у(х) в точ­ках xi , i = 1, 2,..., n.

 

7.2. Метод Эйлера (ломаных)

Рассмотрим уравнение (7.1) в точках xi, i = 0, 1, ... , n - 1 и заменим производную y'(xi) разностной формулой

y'(xi) = (7.4)

Тогда получим рекуррентную формулу метода Эйле­ра для вычисления приближенных значений у(xi + 1):

yi + 1 = yi + h f(хi, уi), i = 0, 1,..., п - 1. (7.5)

Здесь через yi обозначены приближенные значения y(xi), т. е. уi = у (xi), i = 1, 2, ... , п.

На рис. 7.1 дана геометрическая иллюстрация метода Эйлера (7.5). Уравнение касательной к графику решения у(х) в точке (х0, у0) имеет вид

y = y0 + f(x0, y0)(xx0) (7.6)

так как у'(х0) = f(x0, y0). Интегральная кривая у(х) на от­резке [х0, х1] заменяется отрезком касательной (7.6), со­единяющей точку (х0, у0) с точкой (х1, у1), где у1 = у0 + f(x0, у0)(х - х0) (рис. 7.1). Точка (х1, у1) уже не лежит на интегральной кривой у = у(х), удовлетворяющей на­чальному условию (7.2).

При i = 1 формула (7.5) дает точку (х2, у2), которая оп­ределяется с помощью касательной у = y1 + f(x1, у1)(х - x1), проведенной в точке (х1, у1) к интегральной кривой у(х),

 

Рис. 7.1

удовлетворяющей уравнению (7.1) и начальному условию y(x1) = y1.

Таким образом, с каждым шагом i метод Эйлера (7.5) дает точки (xi, yi), которые, вообще говоря, удаляются от интегральной кривой, соответствующей точному решению задачи Коши (7.1), (7.2). Вместо интегральной кривой ме­тод Эйлера дает ломаную, изображенную на рис. 7.1, по­этому его часто называют методом ломаных.

Формулу (7.5) можно получить и другим способом. Рассмотрим разложение искомого решения у(х) в ряд Тейлора в окрестности точки х0:

у(х) = у(х0) + у¢(х0)(х - х0) + + ... .

Ограничившись двумя слагаемыми и учитывая, что y'(х0) = f(x0, у0), при х = х1 получим (7.5).

Также здесь получен следующий результат – погреш­ность вычисления значения у1 есть величина порядка O(h2), а погрешность значения уп – величина порядка O(h). Из-за большой погрешности метод Эйлера применя­ется редко.

 

7.3. Метод Эйлера с уточнением

Для повышения точности метода Эйлера применяют следующий прием. Сначала находят приближенное зна­чение решения по методу Эйлера:

,

а затем уточняют его по формуле

.

Этот метод называется методом Эйлера-Коши, и ре­куррентные соотношения для его реализации могут быть записаны в виде

, (7.7)

i = 0, 1, ..., n - 1

Метод Эйлера-Коши имеет погрешность порядка O(h2). Геометрическая иллюстрация метода Эйлера-Коши пока­зана на рис. 7.2. Очередное приближение метода Эйлера-Коши соответствует точке пересечения диагоналей парал­лелограмма, построенного на двух звеньях ломаной метода Эйлера.

 

Метод Эйлера-Коши является одним из частных слу­чаев методов Рунге-Кутта.

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

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

Выведем формулы метода Рунге-Кутта второго поряд­ка точности. Для этого решение представим куском ряда Тейлора, отбрасывая члены с порядком выше второго. Тогда приближенное значение искомой функции в точке xi, можно записать в виде:

y1 = у(х0) + у¢(х0)(х - х0) + = y0 + f(x0, y0)h + . (7.8)

Вторую производную у"(х0) можно выразить через производную функции f(x, у), однако в методе Рунге-Кут­та вместо производной используют разность

 

,

соответственно подбирая значения параметров , , . Тогда (7.8) можно переписать в виде

y1 = y0 + h[bf(x0, y0) + αf(x0 + γh, y0 + δh)], (7.9)

где α, β, γ и δ – некоторые параметры.

Рассматривая правую часть (7.9) как функцию аргу­мента h, разложим ее по степеням h:

y1 = y0 + (α + β)hf(x0, y0) + αh2fx(x0, y0) + δf(x0, y0)],

и выберем параметры α, β, γ и δ так, чтобы это разложе­ние было близко к (7.8). Отсюда следует, что

α + β = 1, αγ = 0,5, αδ = 0,5 f(х0, y0).

С помощью этих уравнений выразим β, γ и δ через па­раметр α, получим

y1 = y0 + h[(1 – α) f(x0, y0) + αf(x0 + , y0 + f(x0, y0))],

0 < α £ 1. (7.10)

Теперь, если вместо (х0, у0)в (7.10) подставить (х1, у1),получим формулу для вычисления у2– приближенного значения искомой функции в точке х2.

В общем случае метод Рунге-Кутта применяется на произвольном разбиении отрезка [х0, X]на п частей, т. е. с переменным шагом

x0, x1, ..., xn; hi = xi + 1 xi, xn = X. (7.11)

Параметр α выбирают равным 1 или 0,5. Запишем окончательно расчетные формулы метода Рунге-Кутта второго порядка с переменным шагом для α = 1:

yi + 1 = yi + hi f(xi + , yi + f(xi, yi)), (7.12)

i = 0, 1, ..., n – 1.

и α = 0,5:

yi + 1 = yi + [f(xi, yi)+ f(xi + hi, yi + hif(xi, yi))], (7.13)

i = 0, 1, ..., n – 1.

Дадим геометрическую интерпретацию методов Рун­ге-Кутта (7.12), (7.13).

На рис. 7.3 а иллюстрируется формула (7.12). Снача­ла по методу Эйлера вычисляется приближенное значение в точке хi + hi/2. В этой точке определяется касатель­ная к интегральной кривой, параллельно которой через точку (xi, yi)проводится прямая до точки пересечения с прямой х = xi + 1.Ордината точки пересечения yi + 1 при­нимается за приближенное значение искомого реше­ния в точке xi + 1.

Рис. 7.3, б интерпретирует формулу (7.13). Методом Эйлера вычисляется приближенное значение yi + hif(хi, yi) в точке хi + 1.В этой точке тангенс угла наклона касатель­ной к интегральной кривой равен выражению f(xi + hi, yi + hif(хi, yi)).Через точку (xi, уi) проводится прямая, тан­генс угла наклона которой определяется как полусумма угловых коэффициентов касательных, проведенных через точки (xi, yi) и (xi +1, yi + hif(xi, уi)).За приближенное зна­чение искомого решения в точке хi +1 принимается орди­ната точки пересечения этой прямой с прямой х = xi+1.

Отметим, что метод (7.13) совпадает с методом Эйлера-Коши (7.7).

 

Рис. 7.3

Теорема 7.2. Если правая часть f(х, у)уравнения (7.1) непрерывна и ограничена вместе со своими производными до второго порядка включительно, то решение, получен­ное по формулам (7.12) и (7.13), равномерно сходится к решению задачи (7.1) и (7.2) с погрешностью O(max hi2).

Приведем наиболее употребительные формулы метода Рунге-Кутта, формулы четвертого порядка точности:

(7.14)

7.5. Метод Рунге-Кутта с автоматическим выбором шага. Правило Рунге оценки погрешности

Для метода Рунге-Кутта применимо правило Рунге для оценки погрешности. Пусть у(х; h) - приближенное значение решения в точке х, полученное по формулам (7.12), (7.13) или (7.14) с шагом h, а р - порядок точнос­ти соответствующей формулы. Тогда погрешность R(h) значения у(х; h) можно оценить, используя приближен­ное значение у(х; 2h) решения в точке х, полученное с шагом 2h:

R(h) = (7.15)

где р = 2 для формул (7.12) и (7.13) и р = 4 для (7.14).

Уточненное решение запишем в виде

(7.16)

В алгоритмах с автоматическим выбором шага предва­рительно задают погрешность в виде положительного параметра ε, и на каждом этапе вычисления следующего значения yi + 1 подбирают такой шаг h, при котором вы­полняется неравенство

(7.17)

7.6. Задача Коши для системы дифференциальных уравнений

Метод Рунге-Кутта применим и к задаче Коши для системы m дифференциальных уравнений первого порядкa с m неизвестными функциями

x Î (x0, X) (7.18)

y1(x0) = y1,0, y2(x0) = y2,0,…, ym(x0) = ym,0 (7.19)

Приведем для задачи (7.18), (7.19) расчетные форму­лы метода Рунге-Кутта четвертого порядка. Пусть требу­ется найти систему т функций у1(х), у2(х),..., уm(х), удовлетворяющих в интервале (х0, X) дифференциаль­ным уравнениям (7.18), а в точке х0 - начальным усло­виям (7.19). Предположим, что отрезок [х0, X] разбит на N частей:

xi = х0 + ih, h = .

Тогда каждую l-ю функцию уl(х) можно приближенно вычислять в точках хi +1 по формулам Рунге-Кутта

Kl,1 = fl(xi, y1,i, y2,i,… ,ym,i), i = 1, 2, ... , m,

Kl,2 = fl(xi + , y1,i + , y2,i + ,… ,ym,i + ),

i = 1, 2, ... , m;

Kl,3 = fl(xi + , y1,i + , y2,i + ,… ,ym,i + ),

i = 1, 2, ... , m;

Kl,4 = fl(xi + h, y1,i + , y2,i + ,… ,ym,i + ),

i = 1, 2, ... , m;

yl,i+1 = yl,i + ( + + + ), (7.20)

i = 1, 2, ... , m;

 

Здесь через уl,i обозначено приближенное значение функции уl(х) в точке хi.

Необходимо обратить внимание на порядок вычислений по фор­мулам (7.20). На каждом шаге сначала вычисляются ко­эффициенты Kl,i в следующем порядке:

 

и лишь затем приближенные значения функций y1,i+1, y2,i+1,…, ym,i+1.

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

Задача Коши для дифференциального уравнения n-го порядка

y(n) = f(x, у, у', ... , y(n-1)), x Î (х0, X), (7.21)

y(x0) = y0, y'(x0) = y1,0,…,y(п-1)(x0) = yп - 1,0. (7.22)

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

z0 = y, z1 = y',..., zn - 1 = y(n 1). (7.23)

Учитывая (7.23), из уравнения (7.21) получим систе­му дифференциальных уравнений

(7.24)

Начальные условия (7.22) для функций zl, переписыва­ются в виде

z0(x0) = y0, z1(x0) = y1,0,…,zn-1(x0) = yп - 1,0. (7.25)

Запишем для полученной системы метод Рунге-Кутта:

zl,i+1 = zl,i + ( + + + ), (7.26)

i = 0, 1, 2, ... , N; l = 0, 1,…, n – 1.

Для вычисления коэффициентов , , , имеем следующие формулы:

K0,1 = z1,i,

K1,1 = z2,i,

…………

Kn-1,1 = f(xi, z0,i, z1,i,…, zn-1,i),

K0,2 = z1,i + ,

K1,2 = z2,i + ,

…………………

Kn-1,2 = f(xi + , z0,i + , z1,i + ,…, zn-1,i + ),

K0,3 = z1,i + ,

K1,3 = z2,i + ,

…………………

Kn-1,3 = f(xi + , z0,i + , z1,i + ,…, zn-1,i + ),

K0,4 = z1,i + ,

K1,4 = z2,i + ,

Kn-1,4 = f(xi + , z0,i + , z1,i + ,…, zn-1,i + ).

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

Примерами краевых задач для обыкновенных диффе­ренциальных уравнений являются следующие:

1) Найти функцию u(х), которая удовлетворяет на от­резке [а, b] уравнению

u"(х) = -f(х), (7.27)

а на концах отрезка условиям

u(а) = u (b) = 0. (7.28)

 

Задача (7.27), (7.28) имеет следующее физическое со­держание. Между точками х = а и х = b натянута упру­гая струна, находящаяся под действием внешней изгиба­ющей нагрузки. f(x) - величина нагрузки, а u(х) - прогиб струны в безразмерных единицах.

2) Двухточечная краевая задача для линейного диффе­ренциального уравнения второго порядка:

u"(x) + p(x)u'(x) + q(x)u(x) = f(x), (7.29)

α0u(a) + α1u'(a) = A,

β0u(b) + β1u'(b) = B. (7.30)

а) Метод прогонки

Рассмотрим частный случай задачи (7.29), (7.30):

u"(x) + p(x)u'(x) + q(x)u(x) = f(x), (7.31)

u(0) = а, u(Х) = b. (7.32)

Введем сетку xi = ih, i = 1, 2, ... , N; h = X/N. Обозна­чим через ui приближенные значения u(x) в узлах сетки. Рассмотрим уравнение (7.31) во внутренних узлах хi i = 1, 2,..., N - 1и заменим производную второго поряд­ка разностной формулой (5.11)

(7.33)

Тогда из (7.31), (7.32) получим для определения ui си­стему линейных уравнений

(7.34)

Система (7.34) при р(x) ≥ 0 имеет решение. Система (7.34) представляет собой систему линейных уравнений с трехдиагональной матрицей

(7.35)

и для ее решения применим метод прогонки, который фактически является методом исключения неизвестных Гаусса.

1. Прямой ход прогонки. Запишем первое уравнение (7.35) в виде

. (7.36)

Подставив (7.36) во второе уравнение (7.35) и упрос­тив выражения, увидим, что можно для ui получить формулы

(7.37)

Из уравнения (7.35), учитывая (7.37) при i = N - 2, получим

(7.38)

2. Обратный ход прогонки. После вычисления прогоночных коэффициентов можно найти значения решения задачи. Формула (7.38) дает значение

uN - 1:

uN - 1 = βN - 1 (7.39)

 

Остальные значения вычисляем в обратном порядке:

(7.40)

б) Метод стрельбы (пристрелки)

Рассмотрим метод пристрелки на примере решения линейной краевой задачи (7.31), (7.32):

u"(x) + p(x)u'(x) + q(x)u(x) = f(x),

u(0) = а, u(Х) = b.

Если известны частное решение u(х) неоднородного уравнения

u"(x) + p(x)u'(x) + q(x)u(x) = f(x),

и два линейно независи­мых решения и1(х), u2(х) однородного уравнения

u"(x) + p(x)u'(x) + q(x)u(x) = 0,

то общее решение неоднородного диффе­ренциального уравнения

u"(x) + p(x)u'(x) + q(x)u(x) = f(x),

можно за­писать в виде

u(x) + С1u1(х) + C2u2(x).

Постоянные С1, С2 можно определить из краевых ус­ловий (7.32).

В методе пристрелки используется следующий способ.

Сначала находят частное решение u(х) неоднородного уравнения

u"(x) + p(x)u'(x) + q(x)u(x) = f(x),

удовлетворяющее усло­вию u(0) = а, и частное решение u1(х) однородного урав­нения

u"(x) + p(x)u'(x) + q(x)u(x) = 0,

удовлетворяющее условию u(0) = 0.

Затем общее решение u(х) неоднородного уравнения

u"(x) + p(x)u'(x) + q(x)u(x) = f(x),

удовлетворяющее условию u(0) = а, записывают в виде u(х) + Сu1(х). Остается найти постоян­ную С1 из условия u(x) + Cu1(x) = b.

Приведем сеточный аналог метода пристрелки. Пусть краевая задача приведена к системе уравнений (7.34)

(7.41)

Найдем частные решения неоднородной (7.42) и одно­родной (7.43) систем уравнений:

(7.42)

(7.43)

Выбирая произвольные значения для (при этом должно быть ), находим из (7.42) и (7.43) формулы для вычисления частных решений:

(7.44)

(7.45)

Найдем С из условия и запишем решение:

(7.46)

(7.47)

 

Алгоритм метода пристрелки заключается в том, чтобы выбрать шаг h и выполнить последовательно вычисления по формулам (7.44) - (7.47).

<== предыдущая лекция | следующая лекция ==>
Численное интегрирование | Уравнения в частных производных

Дата добавления: 2020-08-31; просмотров: 446;


Поиск по сайту:

Воспользовавшись поиском можно найти нужную информацию на сайте.

Поделитесь с друзьями:

Считаете данную информацию полезной, тогда расскажите друзьям в соц. сетях.
Poznayka.org - Познайка.Орг - 2016-2024 год. Материал предоставляется для ознакомительных и учебных целей.
Генерация страницы за: 0.085 сек.