Элементы вычислительной математики Учебное пособие по курсу Информатика

8. Решение Обыкновенных дифференциальных уравнений.

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

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

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

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

. (8.1)

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

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

, (8.2)

где С – постоянная интегрирования, то говорят, что уравнение проинтегрировано в элементарных функциях.

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

, (8.3)

то говорят, что уравнение проинтегрировано в квадратурах.

Если удается получить решение в виде (8.2) или (8.3), то говорят, что уравнение интегрируемо в конечном виде. Многие практические задачи не позволяют найти окончательные решения в конечном виде, поэтому для их поиска приходится применять численные методы.

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

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

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

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

Пусть рассматривается начальная задача

. (8.4)

Для простоты изложения выберем постоянный шаг изменения параметра x и обозначим значения аргумента

а приближенные значения y(x_k) обозначим y_k. Чтобы найти эти значения, заменим в уравнении производную ее приближенным разностным отношением

, т. е.,

откуда

. (8.5)

По формуле (8.5) можно, начиная от y₀ и полагая последовательно k=0, 1, 2, …, найти значения

Метод Эйлера имеет простой геометрический смысл, показанный на рис 8.2, где изображены также интегральные линии. Он состоит в том, что через заданную точку M₀(x₀,y₀) мы проводим не искомую интегральную линию, которая нам не известна, а отрезок M₀M₁ касательной к этой линии в точке M₀. Через M₁ мы проводим отрезок M₁M₂, угол наклона которого соответствует касательной к интегральной линии в точке M₁, и т. д. Следует заметить, что поскольку уже точка M₁ в общем случае не лежит на исходной интегральной кривой, то направление касательной в точке M₁ также не будет полностью совпадать с направлением касательной к исходной интегральной кривой при x=x₁.

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

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

К недостаткам метода следует отнести малую точность и систематическое накопление ошибок. По этим причинам в таком виде метод Эйлера почти не используется.

8.3. Модификации метода Эйлера.

8.3.1. Усовершенствованный метод ломаных.

Идея метода показана на рис.8.3. Согласно классическому методу Эйлера, для определения значения функции в точке x=x₁ необходимо вычислить производную (тангенс угла наклона касательной) в точке x=x₀ и по этому направлению сделать один шаг. Новое значение функции будет соответствовать точке M₁.

По усовершенствованному методу ломаных движение по касательной в точке x=x₀ сначала производится только на половину шага h/2 до точки M₁₀, в которой вычисляется новое направление касательной, и по этому направлению из точки M₀ делается шаг до точки M₁₁.

Метод требует вычисления промежуточного значения функции в точке M₁₀ по формуле

(8.6)

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

. (8.7)

8.3.2. Усовершенствованный метод Эйлера-Коши.

Идея метода иллюстрируется рис.8.4. По этому методу сначала определяется тангенс угла наклона касательной в точке M₀, равный величине , затем по направлению касательной производится полный шаг до точки, где- «грубое» приближение значения функции на исходной интегральной кривой приx=x₁. В точке M₁ определяется направление касательной и вычисляется среднее значение тангенса угла наклона на шагеh . По этому направлению из точки M₀ делается полный шаг h до точки M₁₁, которая и принимается в качестве точки на исходной интегральной кривой.

Метод требует вычисления первого приближения значения функции в точке x=x₁ по формуле Эйлера

и последующего уточнения по формуле

. (8.8)

Усовершенствованный метод ломаных и усовершенствованный метод Эйлера-Коши имеют порядок ошибки h², при уменьшении величины шага в 10 раз суммарная ошибка сокращается приблизительно в 100 раз.

Итерационная обработка по усовершенствованному методу Эйлера-Коши.

Усовершенствованный метод Эйлера-Коши можно уточнить, применяя итерационную обработку каждого значения y_k. Для этого сначала по формуле (8.5) находят первое приближение значения (8.5) находят первое приближение значения у_k₊₁:

которое затем уточняется по формуле (8.8):

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

(8.9)

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

Усовершенствованный метод Эйлера-Коши с итерационной обработкой обеспечивает порядок ошибки h³, сокращение величины шага в 10 раз в среднем сокращает общую погрешность расчета в 1000 раз.

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

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

В точке (x_n, y_n) вычисляется тангенс угла наклона k₁.
Используя его, идем на половину шага вперед и вычисляем тангенс угла наклона k₂.
Используя k₂, опять двигаемся из точки (x_n, y_n) на половину шага и вычисляем новое значение тангенса угла наклона k₃.
Используя k₃, опять начинаем из точки (x_n, y_n), но делаем теперь полный шаг вперед и там вычисляем новый тангенс угла наклона k₄.
Четыре тангенса угла наклона усредняем с весами и, используя это среднее значение, делаем окончательный шаг от (x_n, y_n) к (x_n₊₁, y_n₊₁).

Этот алгоритм соответствует системе уравнений (8.10).

(8.10)

Блок-схема метода приведена на рис.8.5.

Эта схема обеспечивает порядок ошибки h⁴, уменьшение шага в 10 раз снижает общую ошибку приблизительно в 10000 раз.

Метод Рунге-Кутта легко распространяется на системы дифференциальных уравнений с помощью формальной замены переменных. Например, для системы из двух уравнений

(8.11)

получаем следующую расчетную схему:

(8.12)

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

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

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

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

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

Алгоритм расчета, основанный на комбинации формул Рунге-Кутта, был предложен в 1970 году Е. Фельдбергом. Алгоритм требует шести вычислений функции за шаг:

(8.13)

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

. (8.14)

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

. (8.15)

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

Блок-схема метода приведена на рис.8.6. Коэффициенты увеличения и уменьшения шага, указанные в блок-схеме, выбраны по результатам анализа решения конкретной задачи газодинамики. В зависимости от вида функции они могут принимать другие значения. Для того, чтобы обеспечить более гибкий выбор шага, рекомендуется задавать значения коэффициентов несколько отличными друг от друга.

111

С.И. Шувалов, С.Г. Ушаков, О.Б. Сперанская

Содержание