Г Л А В А   6

 

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

 

 

 

§ 1. Основные понятия

 

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

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

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

Рассмотрим простой пример. Пусть степень изменения величины y по отношению к изменению x пропорциональна y. Математически это утверждение запишется в виде простейшего дифференциального уравнения

 

.                                                (6.1)

 

Решение этого уравнения хорошо известно:

 

,

 

где a – произвольная постоянная. При различных значениях постоянной a получается семейство кривых, которые все удовлетворяют уравнению (6.1). Собственно говоря, уравнение (6.1) является просто утверждением, что в каждой точке кривой значение самой функции равно значению производной. Если в дополнение к дифференциальному уравнению задать значение y для некоторого значения x, то можно определить постоянную a. Например, предположим, что решение уравнения (6.1) должно проходить через точку , , что обычно записывается в следующем виде:

 

.                                              (6.2)

 

При этом легко найти, что постоянная a равна 1 и что из всего семейства кривых только одна удовлетворяет одновременно (6.1) и (6.2):

 

.

 

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

 

,                             (6.3)

 

где x – независимая переменная.

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

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

 

.                                    (6.4)

 

Такая форма записи называется уравнением, разрешенным относительно старшей производной. Примером такого уравнения является рассмотренное ранее уравнение (6.1)

Решением дифференциального уравнения n-го порядка (6.3) называется всякая n раз дифференцируемая функция , которая после ее подстановки в уравнение обращает его в тождество. График функции  называют интегральной кривой.

Общее решение обыкновенного дифференциального уравнения (6.3) содержит n произвольных постоянных :

 

.                                   (6.5)

Функция (6.5) является решением уравнения (6.3) при любых значениях постоянных , т.е. уравнение (6.3) имеет бесконечное множество решений. Единственные (частные) решения получают с помощью дополнительных условий, которым должны удовлетворять искомые решения. При этом постоянные  получают конкретные значения, определяемые дополнительными условиями.

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

Задача Коши. Если дополнительные условия задаются в одной точке, то такая задача называется задачей Коши. Дополнительные условия в задаче Коши называются начальными условиями, а точка , в которой они задаются, – начальной точкой. Для уравнения (6.3) начальные условия могут быть заданы следующим образом:

 

   …, 

 

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

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

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

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

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

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

Наиболее распространенным и универсальным подходом к численному решению дифференциальных уравнений является метод конечных разностей. Суть этого подхода состоит в следующем. Область непрерывного изменения независимой переменной (аргумента) заменяется дискретным множеством точек, называемых узлами. Эти узлы образуют расчетную сетку. Искомая функция непрерывного аргумента приближенно заменяется функцией дискретного аргумента на заданной сетке. Эта функция называется сеточной. Исходное дифференциальное уравнение заменяется разностным соотношением относительно сеточной функции. При этом входящие в уравнения производные заменяются разностными отношениями (см. гл. 4, §1). Такая замена дифференциального уравнения разностным уравнением называется его аппроксимацией  на сетке. Таким образом, решение дифференциального уравнения сводится к отысканию значений сеточной функции в узлах сетки.

 

 

§ 2. Решение задачи Коши

 

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

 

                                 (6.6)

 

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

,

можно переписать в следующем виде:

,

,

где z – новая зависимая переменная, определяемая вторым уравнением. Теперь получается система уравнений относительно y и z. Решение этой системы даст функцию и ее производную. <

Построение численных алгоритмов решения уравнения (6.6) опирается на дискретизацию задачи. Введем в области расчета  дискретный набор точек , , , в которых будет вычисляться приближенное решение (рис. 6.1). Точки  будем называть узлами интегрирования или узлами сетки, расстояние h между узлами – шагом интегрирования или шагом сетки. Совокупность всех узлов  будем называть сеточной областью или просто сеткой узлов.

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

 – совокупность искомых приближенных значений решения задачи (6.6) в узлах сетки;

 – совокупность значений правой части уравнения (6.6) в узлах.

Различные совокупности величин, отнесенных к узлам сетки, называются сеточными функциями[1].


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

 

 ,

 

где  – значение точного решения в узле сетки.

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

Переходим к обсуждению конкретных методов получения приближенного решения задачи (6.6) в узлах сетки. Простейший способ их конструирования опирается на замену производной в левой части уравнения (6.6) в окрестности каждого узла приближенным разностным отношением по формулам численного дифференцирования (см. гл. 4, §1).

Метод Эйлера. Заменяя в (6.6) производную в окрестности каждого i-го узла сетки правым разностным отношением, приходим к методу Эйлера:

 

                  (6.7)

 

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

Замкнутую систему разностных уравнений вместе с дополнительными условиями (начальными или краевыми) называют разностной схемой. Таким образом, (6.7) – это разностная схема Эйлера.<

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

 

,        (6.8)

 

которая непосредственно следует из верхнего соотношения (6.7).

Метод Эйлера имеет очень простую геометрическую интерпретацию. Искомая интегральная кривая  на отрезке  приближается ломаной (рис. 6.2), наклон которой на каждом элементарном участке  определяется наклоном интегральной кривой уравнения (6.7) в точке .

Ø Замечание. К этому же методу можно придти, заменяя производную в уравнении (6.7) левым разностным отношением

 

 

Последовательные значения  в этом случае вычисляются по формуле

.

 

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

 

,   .

 

Такого рода методы,  в которых для вычисления приближенного решения в очередном i-ом узле необходимо дополнительно решать некоторые уравнения (линейные или нелинейные) называются неявными методами. В противоположность этому методы, в которых приближенное решения в очередном i-ом узле явно выражается через предыдущие значения  называются явными методами. При этом, если для вычисления  используется только одно предыдущее значение , то метод называется одношаговым, а если несколько предыдущих значений – многошаговым. Таким образом, метод Эйлера (6.7) является явным одношаговым методом.<

Оставляя вопрос оценки погрешности методов численного интегрирования за рамками нашего изложения отметим, что рассмотренный метод Эйлера обладает первым порядком точности, т.е. .

Чтобы понять, как можно строить методы, обладающие большей точностью обратимся к рис. 6.3. Здесь в пределах отрезка  в увеличенном масштабе изображена интегральная кривая, выходящая из точки с координатами ;    приближенное значение (при ), которое получается по явному методу Эйлера;  – значение, вычисляемое по неявному методу Эйлера (соответствующий отрезок, проведен из точки  с наклоном, равным наклону касательной к интегральной кривой в точке );  – значение, которое соответствует пересечению  с прямой, проведенной из точки  с наклоном, равным наклону касательной к интегральной кривой в середине отрезка . Исходя из приведенного рисунка можно предположить, что точность  больше, нежели  или . Опираясь на эти простые геометрические соображения, сконструируем другие расчетные схемы.

Модифицированный метод Эйлера. Запишем уравнение прямой выходящей из точки  с наклоном, равным наклону интегральной кривой в середине отрезка  (на рис. 6.3. этой прямой соответствует прямая, проходящая через точки   и ):

.

Для точки  получаем

,

откуда немедленно следует разностное соотношение:

.                          (6.9)

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

.                                  (6.9’)

Исключая из (6.9) промежуточные значения  (подстановкой выражения (6.9’)), можно записать этот метод в виде следующей разностной схемы:

         (6.10)

Система уравнений (6.10) представляет собой явную одношаговую разностную схему, аппроксимирующую исходную задачу (6.6) со вторым порядком точности ().

Метод Эйлера с пересчетом. Построение другого способа расчета опирается на следующие геометрические соображения. С помощью метода Эйлера определяется точка , лежащая на прямой  (рис. 6.4). В этой точке снова вычисляется наклон  интегральной кривой, на рисунке этому значению соответствует прямая . Усреднение двух тангенсов дает прямую . Наконец через точку  проводим прямую , параллельную . Точка, в которой прямая   пересечется с ординатой   и будет искомой точкой . Соответствующая разностная схема будет иметь следующий вид:

 

                 (6.11)

где

 

                             (6.11’)

 

Уравнения (6.11)-(6.11’) аппроксимируют исходное дифференциальное уравнение (6.6) со вторым порядком точности.

Последовательные значения , в соответствии с (6.11), вычисляются по формуле

 

.               (6.11’’)

 

Исключая из (6.11’’) , получим явную формулу для

 

.           (6.12)

 

Алгоритм вычислений, описываемый формулой (6.12), часто называется методом Эйлера с пересчетом.

Ø Замечание. Рассмотренные модификации метода Эйлера относятся к группе методов, называемых методами прогноза и коррекции (предиктор-корректор). Суть этих методов состоит в том, что сначала грубо оценивается решение в некоторой точке отрезка интегрирования, а затем оно уточняется с учетом информации о поведении интегральной кривой. В частности, в методе Эйлера с пересчетом сначала вычисляется приближенное значение  по формуле (6.11’), т.е. осуществляется прогнозирование результата, а затем вычисляется уточненное значение  по формуле (6.11’’) с учетом значения .<

Методы Рунге-Кутта. Рассмотренные выше метод Эйлера и его модификации являются частными случаями однопараметрического семейства схем Рунге-Кутта различного порядка точности.  Так, например, модифицированный метод Эйлера (6.10) и метод Эйлера с пересчетом (6.11) представляют собой частный случай схем Рунге-Кутта второго порядка точности, которые определяются следующим разностным соотношением:

 

,                          (6.13)

 

где a – свободный параметр, а  и  – вспомогательные величины, вычисляемые по формулам

 

,     .      (6.13’).

 

Нетрудно видеть, что при  формулы (6.13) и (6.13’) переходят в формулы модифицированного метода Эйлера (6.10), а при  в формулы метода Эйлера с пересчетом (6.11).

Широкое распространение на практике получил метод Рунге-Кутта четвертого порядка точности. Расчетные формулы этого метода имеют следующий вид:

 

                   (6.14)

 

где  – вспомогательные величины.

Метод Рунге-Кутта (6.14) требует существенно большего объема вычислений по сравнению с методом Эйлера и его модификациями, однако это окупается повышенной точностью, что дает возможность проводить вычисления с бóльшим шагом. Другими словами, для получения результатов с одинаковой точностью в методе Эйлера потребуется значительно меньший шаг, чем в методе Рунге-Кутта (6.14).

Представление о многошаговых методах. Многошаговые методы решения задачи Коши характерны тем, что вычисляемое значение решения в текущем узле зависит от данных не только в одном предыдущем узле, но и в ряде предшествующих. В качестве примера построим простейший многошаговый метод. Заменяя в (6.6) производную в окрестности каждого i-го узла сетки центральным разностным отношением (см. гл. 4, §1), получим  следующую разностную схему:

 

                (6.15)

 

Эта система незамкнута так как значение  не определено. Выразим из (6.15) значение :

.                               (6.16)

 

Видно, что искомое значение  зависит от двух предыдущих значений  и . Для того, чтобы начать вычисления по формуле (6.16) при заданном значении  необходимо каким-то образом доопределить значение . Это можно сделать, например, воспользовавшись методом Эйлера: . Формула (6.16) является простейшим многошаговым (двухточечным) методом решения задачи Коши и обладает вторым порядком точности.

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

.

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

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

 

     или    

 

– метод Эйлера (полученный новым способом).

Заменяя подынтегральную функцию полиномом третьей степени можно получить метод Адамса четвертого порядка:

 

.               (6.17)

 

Любопытен вопрос: какой из двух теперь известных нам методов четвертого порядка точности предпочтительней – метод Адамса (6.17) или метод Рунге-Кутта (6.14). При ответе на этот вопрос нужно принимать следующие соображения: метод Адамса требует меньших затрат (арифметических операций) при определении очередного значения , так как при счете по формуле (6.17) нужно лишь один раз вычислять значение функции ,  другие значения –  – к этому моменту уже вычислены (достаточно их сохранять в памяти компьютера), в то время как метод Рунге-Кутта требует в обязательном порядке вычислять четыре вспомогательных значения функции f (см. формулы (6.14)).

С другой стороны, чтобы начать вычисления по формулам Адамса, необходимо помимо заданного значения  как-то определить (например, по тем же формулам Рунге-Кутта) значения , ,  в первых трех узлах интегрирования.

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

Движение частицы описывается уравнением Ньютона:

 

.                                                 (6.18)

 

Это обыкновенное дифференциальное уравнение второго порядка, которое можно заменить эквивалентной системой двух уравнений первого порядка:

 

,   ,                                  (6.19)

 

где v – новая функция, имеющая в данном случае смысл скорости. Уравнения (6.19) необходимо дополнить начальными условиями (постановка задачи Коши)

 

,   .

 

Решение этой задачи хорошо известно: частица  будет совершать гармонические колебания относительно положения равновесия () по закону , где A – амплитуда колебаний, w – частота колебаний, равная ,  – начальная фаза колебаний. Наличие точного аналитического решения этой задачи дает нам возможность оценить точность приближенных численных результатов.

Перейдем к обсуждению численного решения системы уравнений (6.19). Предположим, что решение уравнений (6.19) ищется на интервале времени . Введем на этом интервале дискретную сетку узлов с шагом : {, , }.

Применение метода Эйлера (6.8) для решения системы уравнений (6.19) приводит к простым расчетным формулам:

 

,    ,     .         (6.20)

 

Применение формулы модифицированного метода Эйлера (6.10) сводится к последовательному вычислению значений функций x и v сначала с шагом , с целью определения наклона интегральных кривых в середине отрезка интегрирования, а затем с шагом :

 

,    ,

,   .

 

Исключая из последних формул промежуточные значения , , получим:

 

,     ,

(6.21)

.

 

Можно показать, что применение метода Эйлера с пересчетом (6.12) к данной системе уравнений дает результат, в точности совпадающий с (6.21).

И, наконец, применим для решения системы уравнений (6.19) формулы метода Рунге-Кутта четвертого порядка (6.14):

 

,    ,    ,   

,

,  ,

,  ,

.

 

Ниже в таблице 1 приведены результаты численных расчетов при начальных условиях ,  для различных моментов времени. Масса частицы m и коэффициент упругости k приняты равными единице. Шаг интегрирования  во всех методах равен 0.1.

Таблица 1.

Момент времени t

Точное
значение
x(t)

Метод Эйлера

Модифицированный метод Эйлера

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

0

1

1

1

1

2

-0.41615

-0.45302

-0.41927

-0.41615

4

-0.65364

-0.80974

-0.64892

-0.65365

6

0.96017

1.28642

0.96363

0.96017

8

-0.14550

-0.17751

-0.15880

-0.14549

10

-0.83907

-1.40885

-0.83095

-0.83908

 

Из приведенных результатов видно, что метод Эйлера обладает наибольшей погрешностью, которая, к тому же, со временем сильно растет, т.е. характеризуется склонностью к накапливанию. Модифицированный метод Эйлера дает существенно лучший результат. На всем интервале  относительная погрешность его результатов не превосходит 10%. Метод Рунге-Кутта показывает очень хорошую точность своих результатов с максимальной относительной погрешность не более 0.01%.

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

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

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

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

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

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


§ 3. Решение краевой задачи

 

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

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

.                           (6.18)

 

Краевая задача состоит в отыскании  решения  уравнения (6.18) на отрезке , удовлетворяющего на концах отрезка условиям

 

,   .                                (6.19)

 

Методы решения краевых задач довольно разнообразные – это и точные аналитические методы, и приближенные, и численные.

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

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

Численные методы решения краевой задачи можно разделить на две группы: сведение (редукция) решения краевой задачи к последовательности решений задач Коши и непосредственное применение конечно-разностных методов.

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

Метод стрельбы (“пристрелки”).  Рассмотрим краевую задачу для уравнения второго порядка, разрешенного относительно второй производной:

                         (6.20)

 

Заметим, что решение обсуждаемой задачи имеет простой геометрический смысл. Нужно найти интегральную кривую (“траекторию”), проходящую через точки А, В (рис. 6.4).  Если в (6.20) заменить условие  на , то мы получим задачу Коши (задачу с начальными условиями, заданными в точке ). Считая решение задачи Коши , зависящим от параметра a, будем искать такое его значение , при котором полученная интегральная  кривая  выходит из точки А и попадает в точку В.  Другими словами, нужно подобрать такое значение угла наклона a, чтобы удовлетворялось уравнение

.     (6.21)

 

Следовательно, для нахождения параметра a необходимо решить нелинейное уравнение вида , где .

Необычность ситуации состоит в том, что функция  задана непривычным пока для нас образом – алгоритмически: чтобы найти значение функции  F  при заданном значении аргумента, надо решить задачу Коши

 

                         (6.22)

 

Но с точки зрения численных методов решения нелинейных уравнений (см. главу 2) не важно как задана функция, достаточно уметь вычислять её значения. Если, например, из каких-то соображений (или в результате предварительных расчетов) известно, что искомое решение лежит между двумя интегральными кривыми  и  с начальными наклонами  и , то простейшим методом решения уравнения (6.21) является метод половинного деления. Полагая , решаем задачу Коши при  и, в соответствии с методом половинного деления, отбрасываем один из отрезков:  или , на котором функция  не меняет знак. Процесс поиска решения прекращается, если разность двух последовательно найденных значений a меньше некоторого наперед заданного малого числа. В этом случае полученное последним решение задачи Коши и будет принято за искомое решение краевой задачи.

Для решения уравнения (6.21) можно использовать и другие методы, например метод Ньютона. Его применение состоит в следующем. Пусть  – начальное приближение к , тогда для уточнения корня можно построить итерационный процесс в соответствии с формулой Ньютона (2.15):

 

,                             (6.23)

 

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

 

,                         (6.24)

 

где  – произвольная малая величина.

С учетом (6.24) имеем

 

,        (6.25)

 

Очевидно, что при использовании метода Ньютона уточнения корня уравнения (6.21) для перехода от k-го приближения к -му необходимо два раза решить задачу Коши (6.22): один раз с  и другой – с .

Описанный алгоритм называется методом стрельбы  вполне оправдано, поскольку в нем как бы проводится “пристрелка” по углу наклона интегральной кривой (“траектории”) в начальной точке. Следует отметить, что этот алгоритм хорошо работает в том случае, если решение  не слишком чувствительно к изменению a ; иначе мы можем столкнуться с неустойчивостью решения.  В тех случаях, когда решения дифференциальных уравнений являются быстро растущими функциями, предпочтительней могут оказаться методы, основанные на непосредственной аппроксимации исходной задачи (см., например, [12]).

 

Вопросы для самоконтроля

1.     В чем заключается постановка начальной и краевой задач для обыкновенных дифференциальных уравнений?

2.     Как осуществить переход от уравнения второго порядка к системе уравнений первого порядка?

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

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

5.     Чем отличаются явные и неявных расчетных схемы?

6.     Чем отличаются методы Рунге-Кутта (6.14) и Адамса (6.17)?

 



[1] Фактически сеточная функция – это функция, заданная таблицей значений {xk, yk}.