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

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

Введение:

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

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

Возникает необходимость использовать численные методы, наиболее известным из которых является метод Рунге — Кутты [1]. Что касается Python, то в публикациях по численным методам, например [2,3], данных по применение Рунге — Кутты крайне мало, а по его модификации — методу Рунге-Кутта-Фельберга вообще нет.

В настоящее время, благодаря простому интерфейсу, наибольшее распространение в Python имеет функцию odeint из модуля scipy.integrate. Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие.

Целью настоящей публикации является сравнительный анализ перечисленных средств численного решения систем дифференциальных уравнений с модифицированным автором под Python методом Рунге-Кутта-Фельберга. В публикации так же приведены решения по краевым задачам для систем дифференциальных уравнений (СДУ).

Краткие теоретические и фактические данные по рассматриваемым методам и программным средствам для численного решения СДУ

Для одного дифференциального уравнения n – го порядка, задача Коши состоит в нахождении функции, удовлетворяющей равенству:

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

Перед решением эта задача должна быть переписана в виде следующей СДУ

(1)

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

Модуль имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (задача Коши). Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат odeint(func, y0, t[,args=(), . ]) Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2. yn], а второй – имя независимой переменной.

Функция func должна возвращать список из n значений функций при заданном значении независимого аргумента t. Фактически функция func(y,t) реализует вычисление правых частей системы (1).

Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений при t=t0.

Третий аргумент является массивом моментов времени, в которые вы хотите получить решение задачи. При этом первый элемент этого массива рассматривается как t0.

Функция odeint() возвращает массив размера len(t) x len(y0). Функция odeint() имеет много опций, управляющих ее работой. Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений ei для каждого значения yi по формуле

Они могут быть векторами или скалярами. По умолчанию

Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Она создает объект ОДУ (тип scipy.integrate._ode.ode). Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы. Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей.

Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций. Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши.

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

При численном решении задачи Коши

(2)

(3)

по известному решению в точке t =0 необходимо найти из уравнения (3) решение при других t. При численном решении задачи (2),(3) будем использовать равномерную, для простоты, сетку по переменной t с шагом т > 0.

Приближенное решение задачи (2), (3) в точке обозначим . Метод сходится в точке если при . Метод имеет р-й порядок точности, если , р > 0 при . Простейшая разностная схема для приближенного решения задачи (2),(3) есть

(4)

При имеем явный метод и в этом случае разностная схема аппроксимирует уравнение (2) с первым порядком. Симметричная схема в (4) имеет второй порядок аппроксимации. Эта схема относится к классу неявных — для определения приближенного решения на новом слое необходимо решать нелинейную задачу.

Явные схемы второго и более высокого порядка аппроксимации удобно строить, ориентируясь на метод предиктор-корректор. На этапе предиктора (предсказания) используется явная схема

(5)

а на этапе корректора (уточнения) — схема

В одношаговых методах Рунге—Кутта идеи предиктора-корректора реализуются наиболее полно. Этот метод записывается в общем виде:

(6),

Формула (6) основана на s вычислениях функции f и называется s-стадийной. Если при имеем явный метод Рунге—Кутта. Если при j>1 и то определяется неявно из уравнения:

(7)

О таком методе Рунге—Кутта говорят как о диагонально-неявном. Параметры определяют вариант метода Рунге—Кутта. Используется следующее представление метода (таблица Бутчера)

Одним из наиболее распространенных является явный метод Рунге—Кутта четвертого порядка

(8)

Метод Рунге—Кутта— Фельберга

Привожу значение расчётных коэффициентов метода

(9)

С учётом(9) общее решение имеет вид:

(10)

Это решение обеспечивает пятый порядок точности, остаётся его адаптировать к Python.

Вычислительный эксперимент по определению абсолютной погрешности численного решения нелинейного дифференциального уравнения с использованием обеих функций def odein(),def oden() модуля scipy.integrate и адаптированного к Python методов Рунге—Кутта и Рунге—Кутта— Фельберга

Адаптированные к Python методы Рунге—Кутта и Рунге—Кутта— Фельберга имеют меньшую абсолютную, чем решение с применением функции odeint, но большую, чем с использованием функции edu. Необходимо провести исследование быстродействия.

Численный эксперимент по сравнению быстродействия численного решения СДУ при использовании функции ode с атрибутом dopri5 (метод Рунге – Кутты 5 порядка) и с использованием адаптированного к Python метода Рунге—Кутта— Фельберга

Сравнительный анализ проведём на примере модельной задачи, приведенной в [2]. Чтобы не повторять источник, приведу постановку и решение модельной задачи из [2].

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

где – радиус вектор движущегося тела, – вектор скорости тела, – коэффициент сопротивления, вектор силы веса тела массы m, g – ускорение свободного падения.

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

К системе следует добавить начальные условия: (h начальная высота), . Положим . Тогда соответствующая система ОДУ 1 – го порядка примет вид:

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

Flight time = 1.2316 Distance = 5.9829 Height =1.8542
Flight time = 1.1016 Distance = 4.3830 Height =1.5088
Flight time = 1.0197 Distance = 3.5265 Height =1.2912
Flight time = 0.9068 Distance = 2.5842 Height =1.0240
Время на модельную задачу: 0.454787

Для реализации средствами Python численного решения СДУ без использования специальных модулей, мною была предложена и исследована следующая функция:

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6

Функция increment(f, t, y, tau) обеспечивает пятый порядок численного метода решения. Остальные особенности программы можно посмотреть в следующем листинге:

Время на модельную задачу: 0.259927

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

Решение краевой задачи с поточно разделёнными краевыми условиями

Приведем пример некоторой конкретной краевой задачи с поточно разделенными краевыми условиями:

(11)

Для решения задачи (11) используем следующий алгоритм:

1. Решаем первые три неоднородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

2. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

3. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

4. Общее решение краевой задачи (11) при помощи решений задач Коши записывается в виде линейной комбинации решений:

где p2, p3 — некоторые неизвестные параметры.

5. Для определения параметров p2, p3, используем краевые условия последних двух уравнений (11), то есть условия при x = b. Подставляя, получим систему линейных уравнений относительно неизвестных p2, p3:
(12)
Решая (12), получим соотношения для p2, p3.

По приведенному алгоритму с применением метода Рунге—Кутта—Фельберга получим следующую программу:

y0[0]= 0.0
y1[0]= 1.0
y2[0]= 0.7156448588231397
y3[0]= 1.324566562303714
y0[N-1]= 0.9900000000000007
y1[N-1]= 0.1747719838716767
y2[N-1]= 0.8
y3[N-1]= 0.5000000000000001
Время на модельную задачу: 0.070878

Вывод

Разработанная мною программа отличается от приведенной в [3] меньшей погрешностью, что подтверждает приведенный в начале статьи сравнительный анализ функции odeint с реализованным на Python метода Рунге—Кутта—Фельберга.

3. Н.М. Полякова, Е.В. Ширяева Python 3. Создание графического интерфейса пользователя (на примере решения методом пристрелки краевой задачи для линейных обыкновенных дифференциальных уравнений). Ростов-на-Дону 2017.

Численные методы решения краевых задач

Постановка задачи и основные положения

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

и краевые условия

где [math]F \bigl(x,y,y’,\ldots,y^<(n)>\bigr);

j=\overline[/math] — функции указанных аргументов, заданные в некоторой области их изменения; [math]L[/math] и [math](n-L)[/math] — число условий на левом и правом концах отрезка [math][a,b][/math] соответственно. Общее количество условий равно порядку дифференциального уравнения. Требуется найти функцию [math]y=y(x)[/math] , которая на отрезке [math][a,b][/math] удовлетворяет уравнению (7.1), а на концах отрезка — краевым условиям (7.2).

Если уравнения (7.1),(7.2) линейны относительно искомой функции и ее производных, то краевая задача называется линейной.

Для простоты ограничимся частным случаем линейной краевой задачи для дифференциального уравнения второго порядка [math](n=2)[/math] , которая наиболее часто ставится в вычислительной практике и записывается в виде

(\Omega \equiv [a,b]),[/math]

где [math]p(x),\, q(x),\, f(x)\in C_2[a,b][/math] — заданные функции, а [math]\alpha_0,\,\alpha_1,\, \beta_0,\, \beta_1,\,A,\,B[/math] — заданные числа, 0,

j=0;1[/math] . Требуется найти функцию [math]y(x)[/math] , удовлетворяющую уравнению (7.3) и краевым условиям (7.4). Краевые условия при [math]\alpha_\ne0,

j=0;1[/math] , задают линейную связь между значениями искомого решения и его производной на концах отрезка [math][a,b][/math] .

В простейшем случае, когда [math]\beta_0=0,

\beta_1=0[/math] , краевые условия задают на концах отрезка [math][a,b][/math] только значения функции [math]y(a),\,y(b)[/math] . Такие функциональные условия называют краевыми условиями первого рода. В этом случае краевая задача называется первой краевой задачей.

В случае, когда [math]\alpha_0=0,

\alpha_1=0[/math] , т.е. на концах отрезка заданы только значения производных, краевые условия являются дифференциальными. Такие краевые условия называют условиями второго рода или «мягкими». Последнее название обусловлено тем, что они определяют на концах отрезка [math][a,b][/math] всего лишь наклоны интегральных кривых, а не значения функции [math]y(x)[/math] . В этом случае задача (7.3),(7.4) называется второй краевой задачей.

В общем случае, когда [math]\alpha_0[/math] и (или) [math]\alpha_1;

\beta_0[/math] и (или) [math]\beta_1[/math] не равны нулю, краевые условия носят функционально-дифференциальный характер и называются условиями третьего рода. Тогда задача (7.3),(7.4) называется третьей краевой задачей.

Например, условия [math]y(a)=A,

y(b)=B[/math] являются условиями первого рода. Геометрически это означает, что при решении первой краевой задачи требуется найти интегральную кривую уравнения (7.3), проходящую через данные точки [math](a,A),\, (b,B)[/math] (рис. 7.1,а). Условия [math]y'(a)=A,\, y'(b)=B[/math] являются условиями второго рода. Геометрически вторая краевая задача сводится к отысканию интегральной кривой уравнения, пересекающей прямые [math]x=a,

x=b[/math] под заданными углами [math]\alpha,\,\beta[/math] , где [math]\operatorname\alpha=A,

\operatorname\beta=B[/math] (рис. 7.1,6). Условия [math]y'(a)=A,

y(b)=B[/math] являются частным случаем краевых условий третьего рода, так как [math]\alpha_0=0,

\beta_1=0[/math] . Геометрически данная краевая задача сводится к отысканию интегральной кривой уравнения, проходящей через точку [math](b,B)[/math] и пересекающей прямую [math]x=a[/math] под данным углом [math]\alpha[/math] , где [math]\operatorname\alpha= A[/math] (рис. 7.1,в).

В общем случае краевая задача может:

а) иметь единственное решение;

б) не иметь решений;

в) иметь несколько или бесконечно много решений.

Утверждение 7.1 (о существовании и единственности решения краевой задачи (7.3),(7.4)). Для того чтобы существовало единственное решение краевой задачи (7.3),(7.4), необходимо и достаточно, чтобы однородная краевая задача

имела только тривиальное решение [math]y(x)\equiv0[/math] .

Пример 7.1. Найти аналитическое решение следующих краевых задач:

0 \leqslant x \leqslant \frac<\pi><2>,

y\! \left(\frac<\pi><2>\right)-y’\! \left(\frac<\pi><2>\right)=2[/math] (третья краевая задача);

0 \leqslant x \leqslant 1,

y(1)=0[/math] (первая краевая задача).

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

1. Найдем общее решение однородного уравнения [math]y»+y=0[/math] , одинакового для обеих рассматриваемых задач. Так как характеристическое уравнение [math]\lambda^2+1=0[/math] имеет комплексные сопряженные корни [math]\lambda_<1,2>=\pm i= \alpha\pm \beta i[/math] [math](\alpha=0,

\beta=1)[/math] , то общее решение будет

2. Частные решения неоднородных уравнений находятся методом подбора. Подставляя [math]y_<\text>(x)=C[/math] в уравнение [math]y»+y=1[/math] , а [math]y_<\text>(x)=Dx[/math] в уравнение [math]y»+y=-x[/math] , получаем [math]C=1,

D=-1[/math] . Поэтому [math]y_<\text>(x)=1[/math] в случае «а», [math]y_<\text>(x)=-x[/math] в случае «б».

3. Найдем общее решение неоднородного уравнения как сумму общего решения однородного уравнения и частного решения неоднородного уравнения:

а) [math]y(x)=C_1\cos x+C_2\sin x+1[/math] ; б) [math]y(x)=C_1\cos x+C_2\sin x-x[/math] .

4. Определим значения произвольных постоянных из краевых условий третьего рода (случай «а») и первого рода (случай «б»):

а) найдем [math]y'(x)=-C_1\sin x+C_2\cos x[/math] . Тогда

Отсюда [math]C_1=1[/math] и [math]y(x)=1+\cos x[/math] — решение краевой задачи «а»;

б) общее решение [math]y(x)=C_1\cos x+C_2\sin x-x[/math] и, следовательно, [math]y(0)=C_1=0,

y(1)=C_1\cos1+ C_2\sin1-1=0[/math] , отсюда [math]C_2= \frac<1><\sin1>[/math] и [math]y(x)=\frac<\sin x><\sin1>-x[/math] — решение краевой задачи «б». Таким образом, решение краевой задачи представляет собой такое частное решение, которое удовлетворяет краевым условиям.

Рассмотренный метод нахождения аналитического решения краевых задач применим для ограниченного класса задач. Поэтому в вычислительной практике используются численные и приближенно-аналитические методы, позволяющие найти приближенное решение краевых задач, точные аналитические решения которых не могут быть найдены.

Метод сеток

Рассмотрим линейную краевую задачу с краевыми условиями первого рода (первую краевую задачу):

где [math]p(x),q(x),f(x)\in C_2[a,b][/math] — заданные функции; [math]A,\,B[/math] — заданные числа.

Очевидно, любой отрезок [math][a,b][/math] , на котором ищется решение краевой задачи, можно привести к отрезку [math][0;1][/math] с помощью линейного преобразования [math]\widetilde= \frac[/math] . Действительно, тогда новая переменная [math]\widetilde\in [0;1][/math] . В результате без ограничения общности краевая задача (7.5) может быть решена сначала на отрезке [math][0;1][/math] , а затем это решение с помощью преобразования [math]x=a+(b-a)\cdot \widetilde[/math] может быть записано на отрезке [math][a,b][/math] . То же относится и к исследованию свойств полученного решения.

Утверждение 7.2 (о единственности решения краевой задачи (7.5)). Если функции [math]p(x),q(x),f(x)[/math] принадлежат классу [math]C_2[a,b],

q(x) \geqslant 0[/math] на [math][0;1][/math] , то краевая задача (7.5) имеет единственное решение [math]y(x)\in C_4[0;1][/math] .

Для решения задачи (7.5) применим метод сеток, получаемый путем аппроксимации первой и второй производных. Введем равномерную сетку (где [math]n[/math] — число отрезков разбиения)

Функции [math]p(x),q(x),f(x)[/math] заменяются их проекциями на сетку [math]\Omega_n[/math] , то есть [math]p(x)\to p(x_)=p_i,[/math] [math]q(x)\to q(x_)=q_i,[/math] [math]f(x)\to f(x_)= f_i,[/math] . Вместо точного решения [math]y(x)[/math] отыскивается некоторое приближение [math]\widehat_= \widehat(x_)\approx y(x_),

i=\overline<0,n>[/math] . Первая и вторая производные аппроксимируются на трехточечном шаблоне [math](x_,x_,x_)[/math] по формулам второго порядка (5.10),(5.14):

Краевые условия для этой задачи аппроксимируются точно, т.е. [math]y(a)[/math] и [math]y(b)[/math] заменяются на [math]\widehat_<0>[/math] и [math]\widehat_[/math] . После замены от дифференциальной задачи (7.5) переходим к разностной схеме:

представляющей собой систему алгебраических уравнений трехдиагонального вида:

\delta_=f_[/math] . Здесь система (7.6) записана для внутренних узлов сетки [math]\Omega_n[/math] . Она является трехдиагональной системой линейных алгебраических уравнений и решается методом прогонки.

1. Изложенный метод сеток допускает обобщение. Например, его можно применять для решения нелинейной краевой задачи:

где [math]F(x,y)[/math] — нелинейная по [math]y[/math] функция (в общем случае, который здесь не рассматривается, функция [math]F[/math] зависит также и от [math]y'[/math] ).

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

В силу нелинейности правой части полученная алгебраическая система является нелинейной и для ее решения нельзя использовать метод прогонки в том виде, в каком он изложен для линейной задачи. Поэтому для ее решения используем метод простых итераций, с помощью которого при фиксированном [math]k[/math] (номер итерации) система алгебраических уравнений (7.8) превращается в линейную, так как величины, входящие в правую часть системы, известны из предыдущей итерации. Действительно, для k-й итерации получается система (которая решается на каждой итерации методом прогонки)

Можно показать, что итерации сходятся при выполнении условия [math]q=\frac<1><8>(x_n-x_0)^2M_1 [math]M_1=\max_<[a,b]>\left|\frac<\partial F><\partial y>\right|[/math] с линейной скоростью.

2. Краевые условия второго и третьего рода в задаче, аналогичной (7.5), могут быть аппроксимированы несколькими способами.

Первый способ. Использование аппроксимационных формул (5.4) первого порядка

В силу первого порядка этих аппроксимаций метод сеток в этом случае также будет иметь первый порядок аппроксимации.

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

Третий способ. Применение левосторонней (5.8) и правосторонней (5.9) формул, аппроксимирующих производные со вторым порядком:

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

Алгоритм применения метода сеток

1. Задать сетку [math]\Omega_n[/math] на отрезке [math][a,b][/math] или сформировать ее из условий достижения требуемой точности.

2. Используя аппроксимационные формулы (5.10),(5.14) и один из трех способов аппроксимации краевых условий (в случае, если они второго или третьего рода), перейти от исходной дифференциальной задачи к системе алгебраических уравнений (разностной схеме), неизвестными в которой являются величины, «близкие» к решению краевой задачи в узлах сетки.

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

Пример 7.2. Найти приближенное решение краевой задачи [math]y»+y=1,

0 \leqslant x \leqslant \frac<\pi><2>,[/math] [math]y'(0)=0,[/math] [math]y\! \left(\frac<\pi><2>\right)-y’\! \left(\frac<\pi><2>\right)=2[/math] при [math]n=3[/math] , используя первый способ аппроксимации краевых условий. Записать разностные схемы для второго и третьего способов при произвольном [math]n[/math] .

В поставленной задаче

Для решения задачи воспользуемся методикой.

1. Так как [math]n=3[/math] , то сетка имеет вид [math]\Omega_3=\[/math] , где [math]x_=ih,

y\! \left(\frac<\pi><6>\right)=y_1,[/math] [math]y\! \left(\frac<\pi><3>\right)=y_2,[/math] [math]y\! \left(\frac<\pi><2>\right)=y_3[/math] . Будем искать приближенные значения [math]\widehat_0,\widehat_1, \widehat_2, \widehat_3[/math] . Проекции функций [math]p(x), q(x), f(x)[/math] на сетку имеют вид [math]p_=0,

2. Составим разностную схему. Согласно (7.6), для внутренних узлов сетки получаем

i=1;2[/math] или [math]\widehat_-(2-h^2)\widehat_+ \widehat_=h^2,

Применим первый способ аппроксимации краевых условий. По формуле (5.4) с учетом условия [math]y'(0)=0[/math] на левом конце имеем

На правом конце [math]y\! \left(\frac<\pi><2>\right)=y_3,

y’\! \left(\frac<\pi><2>\right)=y’_3[/math] , и по второй из формул (7.9) [math]\widehat\,’_<3>= \frac<\widehat_<3>-\widehat_<2>>[/math] . Тогда краевое условие [math]y\! \left(\frac<\pi><2>\right)-y’\! \left(\frac<\pi><2>\right)=2[/math] аппроксимируется выражением

В результате получаем разностную схему первого порядка аппроксимации (трехдиагональную систему линейных алгебраических уравнений)

Сравнивая первое уравнение этой системы с рекуррентным соотношением [math]\widehat_= P_\cdot \widehat_+ Q_[/math] метода прогонки, характеризующим обратный ход, получаем [math]P_0=1,

После этого вычисляются все последующие прогоночные коэффициенты по формулам:

Здесь [math]\alpha_,\beta_,\gamma_[/math] соответствуют коэффициентам левой части полученной алгебраической системы, а [math]\delta_[/math] — правой части.

Далее выполняется обратный ход: [math]\widehat_<3>=Q_3,

\widehat_<2>= P_2\widehat_<3>+ Q_2,

\widehat_<1>= P_1\widehat_<2>+ Q_1[/math] .

Результаты решения краевой задачи приведены в табл. 7.1, в которой последний столбец соответствует точному решению [math]y(x)=1+\cos x[/math] , найденному в примере 7.1.

7.1>>\\\hline i& \alpha_& \beta_& \gamma_& \delta_& P_& Q_& \widehat_& y(x) \\\hline 0& 0&-1,\!0000&-1& 0,\!00000& 1,\!00000& 0& 1,\!8648& 2,\!0000\\\hline 1& 1& 1,\!72584& 1& 0,\!27415& 1,\!37771&-0,\!37770& 1,\!8648& 1,\!8666\\\hline 2& 1& 1,\!72584& 1& 0,\!27415& 2,\!87240&-1,\!87242& 1,\!6277& 1,\!5000\\\hline 3& 1& 0,\!47640&-& 1,\!04200&-& 1,\!21853& 1,\!21853& 1,\!0000\\\hline \end[/math]

В силу того, что краевые условия аппроксимированы с первым порядком относительно [math]h[/math] , в данном случае получена разностная схема первого порядка, так как порядок аппроксимации схемы определяется минимальным порядком аппроксимации дифференциального уравнения и краевых условий.

Воспользуемся вторым способом аппроксимации краевых условий для построения разностной схемы второго порядка аппроксимации. Разложим [math]y(x)[/math] в точке [math]x=x_1[/math] относительно точки [math]x_0[/math] по формуле Тейлора:

Выразим из этого соотношения [math]y'(x_0)[/math] и подставим в него вместо [math]y»(x_0)[/math] выражение [math]y»(x_0)=1-y(x_0)=1-y_0[/math] , определяемое исходным дифференциальным уравнением:

Как показывает это соотношение, дифференциальное условие на левой границе аппроксимируется на двухточечном шаблоне [math](x_0,x_1)[/math] со вторым порядком аппроксимации двухточечным алгебраическим уравнением:

Аналогично получается двухточечное алгебраическое уравнение при / [math]i=n-1[/math] и [math]i=n[/math] . Разложение [math]y(x)[/math] в точке [math]x=x_[/math] относительно точки [math]x_n[/math] по формуле Тейлора имеет вид

Выражая отсюда [math]y'(x_n)[/math] с учетом связи [math]y»(x_n)=1-y(x_n)=1-y_n[/math] , следующей из исходного дифференциального уравнения, получаем

Подставим это выражение в граничное условие:

Таким образом, система линейных алгебраических уравнений в окончательном виде записывается следующим образом:

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

Применим третий способ аппроксимации краевых условий для построения разностной схемы второго порядка. Так, для крайней левой точки используется левосторонняя формула (5.8):

Тогда получается трехточечное алгебраическое уравнение:

Аппроксимация производной [math]y’\! \left(\frac<\pi><2>\right)[/math] в крайней правой точке по правосторонней формуле [math]\widehat\,’_= \frac<1> <2h>\bigl(\widehat_-4\widehat_+ 3\widehat_\bigr)[/math] приводит к трехточечному алгебраическому уравнению:

Тогда в этом случае получается следующая система линейных алгебраических уравнений:

Здесь [math]\widehat_<2>[/math] в первом уравнении и [math]\widehat_[/math] в последнем нарушают ее трехдиагональный характер. В этом случае система приводится к трехдиагональному виду путем исключения [math]\widehat_<2>[/math] и [math]\widehat_[/math] из первых двух и последних двух уравнений системы и после этого решается методом прогонки.

Методы минимизации невязки

Описываемые здесь методы относятся к приближенно-аналитическим и могут применяться при решении достаточно широкого класса задач. На основе одного из приближенно-аналитических методов (метода Галеркина) строится метод конечных элементов, излагаемый в разд. 7.5.

Рассмотрим линейную краевую задачу (7.3),(7.4). Ее решение будем искать в виде

где [math]\varphi_0(x), \varphi_1(x), \ldots, \varphi_m(x)[/math] — элементы заданной системы функций; [math]a_1,\ldots,a_m[/math] — неопределенные коэффициенты. Заданная система функций называется базисной, и ее элементы должны удовлетворять условиям:

а) [math]\varphi_(x)\in C_2[a,b],

б) при любом конечном [math]m[/math] функции [math]\varphi_1(x), \ldots, \varphi_m(x)[/math] линейно независимы на отрезке [math][a,b][/math] ;

в) [math]\varphi_0(x)[/math] удовлетворяет краевым условиям (7.4)

г) [math]\varphi_1(x), \ldots, \varphi_m(x)[/math] удовлетворяют условиям

называется невязкой . Она равна разности левой и правой частей уравнения (7.3), образующейся при подстановке [math]\widehat_(x)[/math] вместо [math]y(x)[/math] в дифференциальное уравнение, и характеризует степень отклонения функции [math]\widehat_(x)[/math] от точного решения краевой задачи. Если при некоторых значениях коэффициентов [math]a_1,\ldots,a_m[/math] невязка тождественно равна нулю на отрезке [math][a,b][/math] , а именно

то функция [math]\widehat_(x)[/math] совпадает с точным решением краевой задачи (7.3),(7.4), так как удовлетворяются и уравнение, и краевые условия.

Однако при решении краевых задач, как правило, не удается получить невязку тождественно равной нулю. Поэтому ставится задача: вычислить коэффициенты [math]a_1,\ldots,a_m[/math] таким образом, чтобы невязка в каком-либо смысле стала меньшей. Полученные в результате коэффициенты определяют приближенное решение (7.11).

Выражение для невязки [math]\varepsilon(x; a_1,\ldots, a_m)[/math] с учетом (7.11) удобно записывать в следующей эквивалентной форме:

где [math]L\widehat_\equiv \widehat\,»_(x)+ p(x)\widehat\,’_(x)-q(x) \widehat_(x),

L[/math] — линейный оператор задачи (7.3),(7.4) (выполняются равенства [math]L(y+z)= Ly+Lz,[/math] [math]L(Cy)=C\cdot Ly[/math] для любых [math]y,\,z[/math] и постоянной [math]C[/math] ).

Рассмотрим различные методы, минимизирующие невязку .

А. Метод коллокации. На интервале [math](a,b)[/math] задаются т точек [math]x_1,\ldots, x_n[/math] (точек коллокации) и требуется, чтобы в каждой из них невязка (7.14) обращалась в нуль:

С учетом (7.16) эта система принимает вид

Если полученная система [math]m[/math] линейных уравнений совместна, то из нее определяются коэффициенты [math]a_1,\ldots, a_m[/math] , которые затем подставляются в (7.11).

Б. Метод наименьших квадратов (непрерывный вариант). Неизвестные коэффициенты [math]a_1,\ldots, a_m[/math] должны обеспечивать минимум интеграла от квадрата невязки:

Для решения задачи применяются необходимые условия безусловного экстремум:

Подставляя (7.16) в (7.19), получаем систему [math]m[/math] линейных алгебраических уравнений для нахождения коэффициентов [math]a_1,\ldots, a_m\colon[/math]

В. Метод наименьших квадратов (дискретный вариант). Неизвестные коэффициенты [math]a_1,\ldots,a_m[/math] должны обеспечивать минимум суммы квадратов значений невязки в заданном наборе точек [math]x_1,\ldots,x_n;

n \geqslant m[/math] , то есть [math]x_\in (a,b),

Для решения задачи применяются необходимые условия безусловного экстремума

Отсюда следует система [math]m[/math] линейных уравнений для нахождения коэффициентов [math]a_1,\ldots,a_m[/math] , которая по форме записи совпадает с (7.20), но скалярное произведение определяется по формуле [math]\textstyle<(f,g)= \sum\limits_^ f(x_)g(x_)>[/math] .

Замечание. При [math]n=m[/math] результаты, полученные точечным методом наименьших квадратов и методом коллокации, совпадают. В этом случае точки [math]x_1,\ldots, x_n[/math] являются точками коллокации.

Г. Метод моментов (взвешенных невязок). Неизвестные коэффициенты ах. ат находятся из условия равенства нулю /и моментов невязки:

j=\overline<1,m>[/math] — функции, удовлетворяющие условиям:

б) функции [math]\psi_(x)[/math] являются элементами системы степеней [math]x[/math] или системы тригонометрических функций.

j=\overline<1,m>[/math] называются весовыми, а условие (7.22) является условием ортогональности невязки к весовым функциям.

Д. Метод Галсркина. Он является частным случаем метода моментов, когда в качестве весовых функций используются базисные. Коэффициенты [math]a_1,\ldots,a_m[/math] находятся из условия ортогональности функций базисной системы [math]\varphi_1(x),\ldots, \varphi_(x)[/math] к невязке:

Отсюда следует система [math]m[/math] линейных уравнений для нахождения коэффициентов:

Известно, что при достаточно большом [math]m[/math] условие (7.23) обеспечивает малость невязки в среднем.

Алгоритм применения методов минимизации невязки

1. В выражении (7.11) выбрать систему базисных функций, задать число [math]m[/math] в зависимости от требуемой точности.

2. Найти коэффициенты [math]a_1,\ldots,a_m[/math] путем решения одной из систем алгебраических уравнений (7.18),(7.20),(7.24) в зависимости от выбранного метода.

3. Выписать приближенное решение краевой задачи по формуле (7.11).

Пример 7.3. Найти приближенное решение краевой задачи [math]y»+y=-x,

0 \leqslant x \leqslant 1,[/math] [math]y(0)=0,

y(1)=0[/math] методом коллокации, интегральным методом наименьших квадратов, методом Галеркина

В поставленной задаче

Точное решение найдено в примере 7.1.

Воспользуемся сначала методом коллокации.

1. Зададим [math]m=2[/math] и будем искать решение в виде

где [math]\varphi_0(x)\equiv0[/math] (эта функция удовлетворяет каждому из краевых условий, т.е. [math]\varphi_0(0)=0,

\varphi_0(1)=0[/math] ), функции [math]\varphi_1(x)= x(1-x),

\varphi_2(x)= x^2(1-x)[/math] . Функции [math]\varphi_1(x),\, \varphi_2(x)[/math] линейно независимые, дважды непрерывно дифференцируемые и удовлетворяют условию (7.13). Действительно,

Таким образом, решение краевой задачи ищется в форме

2. Так как [math]m=2[/math] и [math]\varphi_0(x)\equiv 0[/math] , то система (7.18) имеет вид

Выберем узлы коллокации: [math]x_1=1\!\!\not<\phantom<|>>\,4,

Таким образом, имеем линейную систему относительно [math]a_1[/math] и [math]a_2\colon[/math]

3. Приближенное решение задачи: [math]\widehat_2(x)= \frac<217>(42+40x)[/math] .

Решим теперь задачу методом наименьших квадратов (см. непрерывный вариант).

1. Решение краевой задачи ищется в форме [math]\widehat_2(x)= a_1\cdot x(1-x)+ a_2\cdot x^2(1-x)[/math] .

2. Так как [math]f(x)=-x,

\varphi_0(x)\equiv 0[/math] , то система (7.20) имеет вид

Итак, имеем линейную систему относительно [math]a_1[/math] и [math]a_2\colon[/math]

Приближенное решение задачи: [math]\widehat_2(x)=0,\!1875419x(1-x)+ 0,\!1694707x^2(1-x).[/math] .

Решим задачу методом Галеркина.

1. Пусть сначала [math]m=1[/math] . Решение ищется в форме [math]\widehat_1(x)= a_1\cdot x(1-x)[/math] .

2. Тогда система (7.24) преобразуется к виду

Так как [math]\varphi_1(x)= x(1-x),

L\varphi_1(x)= \varphi»_1(x)+ \varphi_1(x)=-2+x(1-x)[/math] , получаем

После вычисления интегралов имеем уравнение [math]-\frac<3><10>\,a_1=-\frac<1><12>[/math] , откуда [math]a_1=\frac<5><18>[/math] .

3. Приближенное решение краевой задачи: [math]\widehat_1(x)=\frac<5><18>\,x(1-x)[/math] . Пусть теперь [math]m=2[/math] .

1. Решение краевой задачи ищется в форме [math]\widehat_2(x)=a_1\cdot x(1-x)+ a_2\cdot x^2(1-x)[/math] .

2. Тогда система (7.24) имеет вид

Вычисляя интегралы, находим

3. Приближенное решение краевой задачи: [math]\widehat_2(x)= x(1-x)\! \left(\frac<71><369>+ \frac<7><41>\,x\right)[/math] .

Сопоставим полученные решения с точным (табл. 7.2).

7.2>>\\\hline x& y_<\text>& y_<\text>& y_<\text>& \text \\\hline 0,\!25& 0,\!045& 0,\!04311& 0,\!0440& 0,\!044014 \\\hline 0,\!50& 0,\!071& 0,\!06807& 0,\!0698& 0,\!069747 \\\hline 0,\!75& 0,\!062& 0,\!05899& 0,\!0600& 0,\!060050 \\\hline \end[/math]

Очевидно, метод Галеркина дал более точный результат.

Пример 7.4. Найти приближенное решение краевой задачи [math]y»+2xy’-2y=2x^2,

0 \leqslant x \leqslant 1,[/math] [math]y'(0)=-2,

y(1)+y'(1)=0[/math] методом Галеркина.

В поставленной задаче

1. Зададим [math]m=2[/math] и подберем функции [math]\varphi_0(x),\, \varphi_1(x),\, \varphi_2(x)[/math] , используя систему [math]1,x,x^2,\ldots[/math] . Функция [math]\varphi_0(x)[/math] должна удовлетворять условиям (7.12):

Пусть [math]\varphi_0(x)=b+cx[/math] , где [math]b,\,c[/math] — неопределенные коэффициенты. Тогда

Отсюда [math]b=4[/math] и [math]\varphi_0(x)=4-2x[/math] .

Функции [math]\varphi_1(x),\, \varphi_2(x)[/math] должны удовлетворять условиям (7.13):

Первое условие выполняется для функций вида [math]\varphi_= x^+b_[/math] . Значения [math]b_[/math] находятся из второго условия [math]1+b_+j+1=0[/math] , откуда [math]b_=-j-2[/math] . Тогда получаем [math]\varphi_1(x)=x^2-3,

Таким образом, решение краевой задачи ищется в форме

2. Тогда система (7.24) имеет вид

3. Приближенное решение краевой задачи [math]\widehat_2(x)= x^2-2x+1[/math] .

Методы сведения краевой задачи к задаче Коши

Метод стрельбы. Суть этого метода заключается в сведении решения краевой задачи к многократному решению задачи Коши. Принцип построения метода стрельбы рассмотрим на примере нелинейной краевой задачи:

где [math]f(x,y,y’)[/math] — нелинейная функция, обусловливающая нелинейность дифференциального уравнения (7.25).

При введении новой переменой [math]z=y'[/math] уравнение (7.25) записывается в нормальной форме Коши, а краевые условия видоизменяются:

где [math]\eta=y'(a)=\operatorname\alpha[/math] — параметр, равный тангенсу угла наклона интегральной кривой в точке [math]x=a[/math] . Угол [math]\alpha[/math] (параметр [math]\eta[/math] ) в процессе многократного решения краевой задачи должен принять такое значение, чтобы интегральная кривая «попала в цель», т.е. в точку [math](b,B)[/math] (рис.7.2 ,а). В общем случае полученное при некотором значении [math]\eta[/math] решение [math]y(x,\eta)[/math] не будет удовлетворять условию [math]y(b,\eta)=B[/math] на правом конце отрезка.

Следовательно, требуется найти такое значение параметра [math]\eta[/math] , чтобы оно было корнем нелинейного уравнения [math]\Phi(\eta)= y(b,n)-B=0[/math] . Для решения этого уравнения, как правило, используются методы половинного деления или секущих. В случае использования метода половинного деления сначала делают «пробные» выстрелы при выбранных наугад или в соответствии с некоторым алгоритмом значениях [math]\eta[/math] до тех пор, пока среди значений [math]\Phi(\eta)[/math] не окажется двух противоположных по знаку. Им соответствует начальный интервал неопределенности, который далее последовательно сокращается путем деления пополам. При применении метода секущих используется формула

где [math]\eta^<(0)>,\,\eta^<(1)>[/math] — начальные значения параметра, [math]k[/math] — номер итерации. Итерации прекращаются при выполнении условия окончания [math]\bigl|\Phi(\eta^<(k)>)\bigr| \leqslant \varepsilon[/math] или [math]\bigl|\eta^<(k+1)>-\eta^<(k)>\bigr| \leqslant \varepsilon[/math] с некоторым положительным [math]\varepsilon[/math] , характеризующим точность решения задачи.

Замечание. Точность решения краевой задачи зависит не только от точности определения параметра [math]\eta[/math] , но также и от точности решения соответствующей задачи Коши. Поэтому одновременно с уточнением параметра [math]\eta[/math] рекомендуется уменьшать шаг при решении задачи Коши, либо выбирать более точный метод.

Рассмотрим применение метода стрельбы для решения линейной краевой задачи (7.3),(7.4):

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

Pers.narod.ru. Обучение. Лекции по численным методам. Численное решение обыкновенных дифференциальных уравнений

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

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

, где x – независимая переменная, — i-ая производная от искомой функции. n — порядок уравнения. Общее решение ОДУ n–го порядка содержит n произвольных постоянных , т.е. общее решение имеет вид .

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

Ясно, что при n=1 можно говорить только о задачи Коши.

Примеры постановки задачи Коши:

Примеры краевых задач:

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

Численные методы решения задачи Коши для ОДУ первого порядка

Постановка задачи. Найти решение ОДУ первого порядка

на отрезке при условии

При нахождении приближенного решения будем считать, что вычисления проводятся с расчетным шагом , расчетными узлами служат точки промежутка [x₀, x_n].

Целью является построение таблицы

т.е. ищутся приближенные значения y в узлах сетки.

Интегрируя уравнение на отрезке , получим

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

то получим явную формулу Эйлера:

, .

Зная , находим , затем т.д.

Геометрическая интерпретация метода Эйлера:

Пользуясь тем, что в точке x₀ известно решение y(x₀) = y₀ и значение его производной , можно записать уравнение касательной к графику искомой функции в точке : . При достаточно малом шаге h ордината этой касательной, полученная подстановкой в правую часть значения , должна мало отличаться от ординаты y(x₁) решения y(x) задачи Коши. Следовательно, точка пересечения касательной с прямой x = x₁ может быть приближенно принята за новую начальную точку. Через эту точку снова проведем прямую , которая приближенно отражает поведение касательной к в точке . Подставляя сюда (т.е. пересечение с прямой x = x₂), получим приближенное значение y(x) в точке x₂: и т.д. В итоге для i–й точки получим формулу Эйлера.

Явный метод Эйлера имеет первый порядок точности или аппроксимации.

Если использовать формулу правых прямоугольников: , то придем к методу

, .

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

Неявный метод Эйлера имеет первый порядок точности или аппроксимации.

Модифицированный метод Эйлера: в данном методе вычисление состоит из двух этапов:

Данная схема называется еще методом предиктор – корректор (предсказывающее – исправляющее). На первом этапе приближенное значение предсказывается с невысокой точностью (h), а на втором этапе это предсказание исправляется, так что результирующее значение имеет второй порядок точности.

Методы Рунге – Кутта: идея построения явных методов Рунге–Кутты p–го порядка заключается в получении приближений к значениям y(x_i₊₁) по формуле вида