Выбор шага при решении дифференциальных уравнений

Постановка задачи. Типы задач для ОДУ.

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

u (p) (x) = f(x, u, u / , u // , … u (p+1) ) заменой u (k) (x) = y_k(x) можно свести к эквивалентной системе из p уравнений первого порядка, представленных в каноническом виде:

Покажем такое преобразование на примере уравнения Бесселя: .

Предполагая тождественную замену y₁(x Задача Коши. Задается координата u_k(x₀) = u_k0 начальной точки интегральной кривой в (p+1)–мерном пространстве (k = 1…p). Решение при этом требуется найти на некотором отрезке x₀ ≤ x ≤ x_max.

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

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

Рассмотрим задачу Коши для уравнения первого порядка:

В окрестности точки x₀ функцию u(x) разложим в ряд Тейлора:

(4)

Идея этого и последующих методов основывается на том, что если f (x , u) имеет q непрерывных производных, то в разложении можно удержать члены вплоть до O(h q+1 ), при этом стоящие в правой части производные можно найти, дифференцируя (3) требуемое число раз. В случае метода Эйлера ограничимся только двумя членами разложения.

Пусть h – малое приращение аргумента. Тогда (4) превратится в

.

Так как в соответствии с (3) u / (x₀) = f(x₀, u₀), то .

Теперь приближенное решение в точке x₁ = x₀+h можно вновь рассматривать как начальное условие, т.е. организуется расчет по следующей рекуррентной формуле:

(5) ,

где y₀ = u₀, а все y_k – приближенные значения искомой функции (см. рисунок). В методе Эйлера происходит движение не по интегральной кривой, а по касательной к ней. На каждом шаге касательная находится уже для новой интегральной кривой (что и дало название методу – метод ломаных), таким образом ошибка будет возрастать с отдалением x от x₀.

При приближенное решение сходится к точному равномерно с первым порядком точности. То есть, метод дает весьма низкую точность вычислений: погрешность на элементарном шаге h составляет ½ h 2 y // (½(x_k+x_k+1)), а для всей интегральной кривой порядка h 1 . При h = const для оценки апостериорной погрешности может быть применена первая формула Рунге, хотя для работы метода обеспечивать равномерность шага в принципе не требуется.

Метод Эйлера легко обобщается для систем ОДУ. При этом общая схема процесса (5) может быть записана так:

(6),

где i = 1… m – число уравнений, k – номер предыдущей вычисленной точки.

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

Данный метод базируется на предыдущем, однако здесь апостериорная погрешность контролируется на каждом шаге вычисления. Как и в предыдущем случае, рассматриваем задачу Коши на сетке с постоянным шагом h. Грубое значение вычисляется по формуле (5) и затем итерационным циклом уточняется по формуле

, (7)

где m – номер итерации. Итерационный цикл повторяется до тех пор, пока . Данная формула также легко обобщается на решение систем ОДУ. Априорная погрешность метода при m = 1 на каждом шаге порядка h 3 .

Метод Рунге-Кутты II порядка.

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

. (8)

Переходя к приближенному решению y ≈ u и заменяя производные в (8) конечными разностями, получаем в итоге следующее выражение:

, (9)

где 0 ≤ α ≤ 1– свободный параметр. Можно показать, что если f(x,u) непрерывна и ограничена вместе со своими вторыми производными, то решение, полученное по данной схеме, равномерно сходится к точному решению с погрешностью порядка h 2 .

Для параметра α наиболее часто используют следующие значения:

1) α = 1. В этом случае
. Графически это уточнение можно интерпретировать так: сначала по схеме ломаных делается шаг h /2, и находится значение . В найденной точке определяется наклон касательной к интегральной кривой, который и будет определять приращение функции для целого шага, т.е. отрезок [ AB ] (см. рисунок) будет параллелен касательной, проведенной в точке (x_k + h/2, y(x_k+h/2)) к интегральной кривой.

2) α = ½. В этом случае

. Можно представить, что в этом случае по методу Эйлера сначала вычисляется значение функции и наклон касательной к интегральной кривой в точке x_{k +1}. Затем находится среднее положение касательной из сравнения соответствующих наклонов в точках x _k и x_{k +1}, которое и будет использоваться для расчета точки y_k+1.

Метод Рунге-Кутты IV порядка.

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

, где i = 1… p, p – число уравнений в системе.

; k – номер точки, для которой осуществляется расчет;

;

;

.

К достоинствам метода следует отнести высокую точность вычислений. Схемы более высокого порядка точности практически не употребляются в силу своей громоздкости. Также немаловажно, что метод является явным, т.е. значение y_{k +1} вычисляется по ранее найденным значениям за известное заранее число действий.

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

Метод Рунге – Кутта — Мерсона.

Этот метод отличается от метода Рунге – Кутта четвертого порядка возможностью оценивать погрешность на каждом шаге и в зависимости от этого принимать решение об изменении шага. Один из вариантов формул:

;

R_{n +1} = 0.2 k ₄ – 0.3 k ₃ – 0.1 k ₅ — погрешность на каждом шаге.

Пусть задана максимальна погрешность . Если , h = h /2 , и n +1 цикл расчета повторяется (с точки x_n , y_n ) c новым шагом. Если же

, h = 2 h

Автоматический выбор шага позволяет значительно сократить время решения ОДУ.

Схема РКМ обобщается на системы ОДУ аналогично классической схеме Рунге – Кутта.

Метод основан на аппроксимации интерполяционными полиномами правых частей ОДУ.

Пусть с помощью любого из методов, рассмотренных выше, вычислено решение заданного дифференциального уравнения в точках x₁, x₂, x₃ (а в точке x₀ решение и так известно – поставлена задача Коши). Полученные значения функции обозначим как y₀, y₁, y₂, y₃, а значения правой части дифференциального уравнения как f ₀, f₁, f₂, f ₃, где f _k = f (x_k, y_k). Начиная с четвертой точки, на каждом шаге интегрирования дифференциального уравнения вычисления осуществляются по схеме

где P – прогноз решения; Е – вычисление f (x , y ); С – коррекция решения; m – количество итераций коррекции. Схемы такого типа называют «прогноз-коррекция»: это подразумевает сначала приблизительное вычисление решение по формуле низкого порядка, а затем уточнение с учетом полученной информации о поведении интегральной кривой.

Прогноз осуществляется по экстраполяционной формуле Адамса:

. (10)

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

. (11)

Вычисление осуществляется по формуле:

.

Количество итераций m ≤ p, где p – порядок используемого метода. В ходе каждой итерации решается нелинейное уравнение (11) относительно неизвестной y₄ (обычно методом простых итераций).

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

Постановка краевой задачи. Метод стрельбы.

Краевая задача – задача отыскания частного решения системы

(12) , 1 ≤ k ≤ p, на отрезке a ≤ x ≤ b, на котором дополнительные условия налагаются на значения функции u_k(x) более чем в одной точке этого отрезка. В качестве примера моно привести задачу нахождения статического прогиба u(x) нагруженной струны с закрепленными концами:

, a ≤ x ≤ b, u(a) = u(b) = 0. Здесь f (x) имеет смысл внешней изгибающей нагрузки на единицу длины струны, деленной на упругость струны.

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

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

(13)

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

(14) .

Сначала выберем некоторое значение , так что . Это уравнение мы можем решить и определить . Таким образом у нас появились два числа и β, которые будут определять задачу Коши , для системы (13), удовлетворять левому краевому условию, но не удовлетворять второму уравнению (14). Тем не менее, решим систему ОДУ с задачей Коши каким-либо из известных нам численных методов. Получим решение вида

(15) , зависящее от как от параметра.

Теперь мы должны каким-либо способом менять параметр до тех пор, пока не подберется такое значение, для которого будет выполнено условие (16) , т.е. правое краевое условие. Для этого случайным образом берут значения до тех пор, пока среди величин не окажется разных по знаку.

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

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

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

, (17)

где U – температура, – теплоемкость, – коэффициент теплопроводности, q – плотность источников тепла.

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

Рассмотрим в качестве примера одномерную задачу, близкую по смыслу к (17):

. (18)

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

Введем равномерную прямоугольную сетку по x и t с шагом h и τ соответственно и заменить производные соответствующими разностными отношениями. Тогда

. (19)

Здесь 1 ≤ k ≤ N-1 – число точек по координате x; 0 ≤ m ≤ M – число точек по координате t. Число неизвестных в (19) больше числа уравнений, недостающие уравнения выводятся из начальных и граничных условий:

; 0 ≤ k ≤ N.

; ; 0 ≤ m ≤ M.

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

, где 1 ≤ k ≤ N-1.

; . (20)

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

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

Видео:Решение системы дифференциальных уравнений методом ЭйлераСкачать

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

Видео:Метод ЭйлераСкачать

1. Задачи с начальными условиями

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

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

где (f) — известная функция (т.е. нечто, что мы можем вычислить, зная (textbf) и (t) ), (textbf) — состояние системы, (dot<textbf>) — производная (textbf) по времени. Как правило, (textbf) и (dot<textbf>) являются векторами. Как подсказывает само название, в задачах с начальными условиями мы располагаем значением (textbf) в начальный момент времени (t_0) : (textbf(t_0) = textbf_0) . Пользуясь им нужно определить состояние системы (textbf) во все последующие моменты времени.

Рис. 1. Производная функция (f(textbf,t)) определяет векторное поле

Задачу с начальными условиями легко представить наглядно. В двумерном случае (textbf(t)) представляет собой траекторию движения точки (textbf

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

) , если она будет проходить через (textbf) (она может проходить через эту точку, а может и не проходить). Представим, что (f) переносит (textbf

) из точки в точку, подобно океанскому течению. Как только мы поместим (textbf

) на плоскость, «течение» подхватит ее. Место, куда в конечном итоге попадет (textbf

) , зависит от того, куда мы поместили ее в начале, но как только это сделано, все последующее движение точки (textbf

) определяется функцией (f) . Траектория движения (textbf

) под действием (f) называется интегральной кривой векторного поля (см. рис. 2).

Рис. 2. Задача с начальным условием. Стартуя из точки (textbf_0) , частица (textbf

) движется со скоростью, определяемой векторным полем

Мы записали функцию (f) как зависящую от (textbf) и от (t) , но производная функция может как зависеть, так и не зависеть от времени. Если она зависит от времени, то в нашем наглядном представлении движется не только (textbf

) , движется все векторное поле, так что скорость (textbf

) зависит не только от того, где она находится, но и от того, каким путем она попала в это место. Зависимость (dot<textbf>) в этом случае проявляется двойственно: во-первых, вектор производной изменяется сам, во-вторых, точка (textbf

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

Видео:7. Линейные дифференциальные уравнения первого порядка. Метод Бернулли.Скачать

2. Численные решения

Обычно во вводных курсах по теории дифференциальных уравнений рассматривают уравнения, имеющие аналитические решения, т.е. предполагается, что искомую функцию можно записать в виде формулы. Например, для дифференциального уравнения (dot=-kx) , где (dot) означает производную от (x) времени, решением является (x=e^) .

Мы, напротив сосредоточимся исключительно на численных решениях, которые заключаются в выполнении дискретных шагов по времени, начиная с начального момента (textbf(t_0)) . Шаг по времени выполняется так: вычисляется функция (f) , чтобы найти приблизительную величину приращения (textbf) за время (Delta t) — (Deltatextbf) , затем (Deltatextbf) прибавляется к (textbf) , чтобы получить новое состояние системы. Теперь нам известно состояние системы (textbf(t_1)) в момент (t_1 = t+Delta t) . Принимая его за начальное, выполним следующий шаг по времени и найдем состояние системы в момент (t_2 = t+2Delta t) и т.д. Функцию (f) можно рассматривать как черный ящик, на вход которого подаются значения (textbf) и (t) , а на выходе получается значение (dot<textbf>) . Различные методы численного решения используют одно или несколько вычислений производной функции на каждом шаге по времени.

2.1 Метод Эйлера

Простейший метод численного решения ОДУ называется методом Эйлера. Пусть начальное условие в момент времени (t_0) обозначается как (textbf_0=textbf(t_0)) , а оценка (textbf) в момент времени (t_0+h) — как (textbf(t_0+h)) , где (h) — длина шага. Сделаем шаг в направлении, указанном вектором производной, и вычислим (textbf(t_0+h))

— в этом и состоит метод Эйлера.

Чтобы понять как работает метод Эйлера, представьте картинку двумерного векторного поля (рис. 3). Вместо реальной интегральной кривой, точка (textbf

) будет двигаться по ломаной, каждое звено которой соответствует шагу численного решения. Направление звена определяется оценкой вектора (f) , сделанной в начале шага, длина звена равна длине этого вектора, умноженной на (h) .

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

Метод Эйлера прост, однако недостаточно точен. Рассмотрим двумерную функцию (f) , интегральные кривые которой представляют собой концентрические окружности (рис. 4). Тогда точка (textbf

) , движимая (f) , должна все время оставаться на одной и той же окружности. Вместо этого, с каждым шагом расчета по методу Эйлера, (textbf

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

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

) на окружность большего радиуса. Уменьшение шага не решает проблему, а только уменьшает скорость накопления ошибки.

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

Кроме того, метод Эйлера может быть неустойчивым. Рассмотрим функцию (f=-kx) , которая, как мы видели, дает интегральные кривые экспоненциально убывающие к нулю. Для достаточно малых длин шага мы получим приемлемые приближения, однако, если (h > 1/k) , то мы получим (|Delta x| > |x|) , и приближенное решение будет колебаться в окрестности нуля. Если же длина шага (h = 2/k) , то эти колебания будут возрастающими. Как говорят в таких случаях, система «взрывается».

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

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

Как видно, формула метода Эйлера получается, если обрезать этот ряд, сохранив только первые два слагаемых в правой части равенства. Это означает, что метод Эйлера будет точен только в случае, если все производные порядков выше первого будут равны нулю, то есть когда функция (textbf(t)) линейна. Погрешность метода Эйлера равна «обрезанной» части разложения в ряд Тейлора и определяется, в первую очередь, слагаемым ((/) ddot<textbf>(t_0)) . Следовательно, мы можем записать погрешность шага метода Эйлера как (O(h^2)) (читается: «порядка (h) в квадрате»).

Предположим, что мы уменьшили шаг вдвое и теперь он равен (h/2) . Хотя при этом мы получаем ошибку, равную лишь четверти той, что была при длине шага (h) , но теперь нам нужно сделать вдвое больше шагов, чтобы покрыть такой же по времени интервал. Это означает, что погрешность, накапливаемая на интервале от (t_0) к (t_1) , зависит от длины шага (h) линейно — (O(h)) — меньше делается шаг и пропорционально уменьшается погрешность. На практике это приводит к тому, что для достижения хорошей точности нужно сделать слишком много шагов.

2.2 Метод средней точки 2

Если бы мы, наряду с (dot<textbf>) , могли вычислить и (ddot<textbf>) , то можно было бы достичь точности (O(h^3)) вместо (O(h^2)) , сохранив в разложении в ряд Тейлора еще одно слагаемое

Вспомним, что производная по времени (dot<textbf>) задается функцией (f(textbf(t),t)) . Для простоты будем считать, что (f) не зависит явно от времени, так что (dot<textbf> = f(textbf(t))) . Правило дифференцирования сложной функции дает в этом случае

Чтобы избежать вычисления (f^prime) , которое может быть сложным и требовать большого расхода времени, найдем приближенное выражение для (ddot<textbf>) через (f) (не содержащее (f^prime) ) и подставим его в eqref. При этом порядок точности (O(h^3)) сохраниться.

Для этого разложим (f) в ряд Тейлора

Введем в него (ddot<textbf>) , подобрав величину шага равной

так что eqref будет иметь вид

Отсюда, умножая обе части равенства на (h) (что превращает O(h^2) в O(h^3)) и перегруппировывая слагаемые, получим

Подстановка правой части этого равенства в eqref дает формулу для численного решения

В соответствии с нею сначала выполняется шаг по методу Эйлера, а затем вычисляется оценка второй производной в средней точке шага для получения нового значения (textbf) . Поэтому этот метод называют методом средней точки. Погрешность шага этого метода оценивается как (O(h^3)) , но для этого требуется два вычисления производной. Иллюстрация работы метода приведена на рис. 5.

Рис. 5. Метод средней точки относится к методам 2-го порядка, т.е. его погрешность на интервале расчета — (O(h^2)) . Шаг метода состоит из a) выполнения шага по методу Эйлера; b) оценки производной в средней точке шага; c) вычисления новой точки интегральной кривой.

Достигнутая погрешность (O(h^3)) — далеко не предел. Вычислив (f) еще несколько раз во время шага, мы сможем исключить производные более высоких порядков. Самый популярный метод, в котором используется эта идея, называется методом Рунге-Кутта 4-го порядка и дает погрешность на шаге (O(h^5)) (метод средней точки можно назвать методом Рунге-Кутта 2-го порядка). Мы не будем выводить, а просто приведем здесь формулы метода Рунге-Кутта 4-го порядка

Видео:Видеоурок "Нахождение частных решений по виду правой части"Скачать

3. Адаптивный выбор шага

Какой бы метод численного решения мы не использовали, остается вопрос выбора длины шага. В идеале, хотелось бы сделать (h) настолько большим, насколько это возможно, чтобы обеспечить заданную точность расчетов и избежать неустойчивости метода. С фиксированной длиной шага мы можем двигаться вперед не быстрее, чем нам позволит это сделать наиболее трудная для вычислений часть (textbf(t)) . Хотелось бы, чтобы длина шага на интервале расчетов уменьшалась, когда это нужно для обеспечения точности, а затем, когда трудное место пройдено, вновь увеличивалась. Т.е. длина шага должна выбираться адаптивно.

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

Допустим, мы вычислили две оценки величины (textbf(t_0+h)) . Первая ( (textbf_a) ) получена вычислением шага (h) от (t_0) к (t_0 + h) , вторая ( (textbf_b) ) — вычислением двух последовательных шагов длины (h/2) от (t_0) к (t_0 + h) . Обе оценки (textbf_a) и (textbf_b) отличаются от истинного значения (textbf(t_0+h)) не более чем на (O(h^2)) . Следовательно и друг от друга (textbf_a) и (textbf_b) отличаются не более чем на (O(h^2)) .

Определим меру текущей погрешности (e) как

Предположим, что требуется обеспечить погрешность на шаге расчета не более (10^) , а текущая погрешность (e) составляет всего (10^) . Т.к. погрешность пропорциональна (h^2) , мы можем увеличить длину шага до

Напротив, если текущая погрешность составляет (10^) мы должны уменьшить шаг до

Адаптивный выбор шага — очень полезная техника.

Видео:18+ Математика без Ху!ни. Дифференциальные уравнения.Скачать

4. Реализация

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

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

С точки зрения решателя, система с которой он работает представляет собой функцию-«черный ящик» (f(textbf,t)) . Решатель должен иметь возможность вычислить (f) при любых значениях (textbf) и (t) и вернуть обновленные значения (textbf) и (t) после окончания расчета. Чтобы обеспечить выполнение этих операций, объект, представляющий ОДУ, должен уметь обрабатывать следующие запросы решателя:

• Возвратить значение (dim(textbf)) . Поскольку (textbf) и (dot<textbf>) — векторы, решатель должен знать их длину чтобы выделить нужное место в памяти, выполнить векторные операции и т. п.
• Получить/задать значения (textbf) и (t) . Решатель должен возвращать окончательные значения этих величин. К тому же, при использовании многошаговых методов, решатель должен задавать промежуточные значения (textbf) и (t) для вычисления производной.
• Вычислить (f) при заданных (textbf) и (t) .

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

Видео:Задача Коши ➜ Частное решение линейного однородного дифференциального уравненияСкачать

Ссылки

1. W.H. Press, B.P. Flannery, S. A. Teukolsky, and W. T. Vetterling. Numerical Recipes in C. Cambridge University Press, Cambridge, England, 1988.

Другое название задач с начальными условиями — задачи Коши. ↩

Функции, реализующие выбор шага и тот или иной численный метод решения ОДУ. ↩

Например, решатель с адаптивным шагом может вызывать решатель с фиксированным шагом. ↩

Видео:Линейное неоднородное дифференциальное уравнение второго порядка с постоянными коэффициентамиСкачать

Моделирование динамических систем: численные методы решения ОДУ

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

Камень бросают вертикально, без начальной скорости с высоты h = 100 м. Пренебрегая сопротивлением воздуха, определить закон движения камня, как функцию высоты камня над поверхностью Земли от времени. Ускорение свободного падения принять равным 10 м/с 2

Видео:Решение физических задач с помощью дифференциальных уравненийСкачать

1. Формализация задачи

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

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

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

где g — ускорение свободного падения у поверхности Земли. Теперь пришло время составить уравнения движения камня. Помните эти уравнения

Левая их часть нас пока не интересует, а вот в правой стоят суммы проекций сил, приложенных к точке на оси координат. Пусть оси x и y располагаются на поверхности, а ось z направлена вверх перпендикулярно ей. Сила одна единственная, её проекции на оси x и y равны нулю, а на ось z проекция отрицательна, так как сила направлена против направления оси, то есть

Масса камня, очевидно не равна нулю, значит можно спокойно поделить обе части получившихся уравнений на неё

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

То что у нас получилось не много не мало — математическая модель процесса происходящего в задаче. Пафосно, да?

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

Аналитическое решение получить просто, даже не буду заморачиваться, оно такое

А вот как решить это численно? И что это вообще такое — «численно»?

Видео:5 Численное решение дифференциальных уравнений Part 1Скачать

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

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

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

проекция ускорения на ось z равна первой производной проекции скорости на ту же оcь z. Ну классно, тогда

вот вам и уравнение первого порядка. Не всегда этот номер проходит (не буду я сейчас про форму Коши!), но в данном случае всё в порядке. Будем искать не координату а скорость точки. Что дальше-то? А дальше

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

Что получается? А вот что

Мы получили приращение скорости. Отрицательное приращение. Как это так, камень, падая вниз будет разгонятся же! Да, будет. Его скорость, вектор его скорости, будет направлен вниз. А значит проекция этого вектора на ось z будет отрицательной. Всё правильно, мы получаем растущую по абсолютному значению проекцию вектора, направленного вниз. Мы знаем, начальное значение скорости — ноль, а значит

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

а ещё через 0.1 секунды

и ещё через 0.1 секунды

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

То есть, воспользовавшись формулой

мы получили зависимость скорости точки от времени. А всего-то нужно взять значение скорости в текущий момент времени, и добавить к нему приращение, которое скорость получит в новый момент времени, отстоящий от текущего на секунд. Приращение времени называется здесь шагом интегрирования. А приращение мы вычисляем как значение производной искомой функции в текущий момент времени умноженное на шаг. Просто? Да просто конечно. И та формула, которую я написал, имеет название название — явный метод Эйлера для численного решения дифференциальных уравнений. Это так называемая рекуррентная формула, когда новое значение вычисляемой величины зависит от её предыдущего значения.

А что же с высотой точки над Землей? Да всё аналогично, смотрите

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

и по этой формуле добавим в нашу таблицу ещё одну колонку

Хм, ну, во-первых, заметно, что высота меняется у нас уже неравномерно, так как скорость со временем меняется. Теперь наша производная сама зависит от времени. Но уже на первом шаге, мы замечаем неладное — скорость уже есть, а вот высота по прежнему 100 метров. Как так?

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

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

Во первых, можно уменьшить шаг. Скажем в 10 раз, пусть секунды

Уже лучше, погрешность в конце счета не превышает 0,05 метров, и это в 10 раз меньше предыдущего значения. Можно предположить, что уменьшив шаг ещё в 10 раз мы получим ещё более точное решение. Я схитрил, выводя значения только для 10 точек с шагом 0.1, на самом деле, чтобы получить такую таблицу нужны уже 100 итерации а не 10. При шаге 0.001 потребуется уже тысяча итераций, а результат будет таким

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

Метод Эйлера самый простой из известных методов интегрирования дифференциальных уравнений. Из нашего простого примера видно, что погрешность метода прямо пропорциональна шагу интегрирования, и это действительно так. Такие методы называются методами 1-го порядка точности.

Точность расчетов даже на шаге 0.1 можно улучшить, если мы применим, скажем метод Рунге-Кутты 4-го порядка точности. Но это отдельная история.

Видео:Линейное дифференциальное уравнение Коши-ЭйлераСкачать

Заключение

Задумайтесь… Мы рассмотрели очень простой пример. Мы даже не применяли компьютер, но уже понимаем принцип, по которому работают те самые мощные в мире суперкомпьютеры, что моделируют ранние этапы жизни Вселенной. Конечно, там всё устроено гораздо сложнее, но принцип лежит этот же самый.

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

📽️ Видео

5 Численное решение дифференциальных уравнений Part 1Скачать

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

МЗЭ 2022 Численное решение дифференциальных уравнений. Неявный метод Эйлера. Ложкин С.А.Скачать

6. Особые решения ДУ первого порядкаСкачать

Дифференциальные уравнения. 11 класс.Скачать

Откуда появляются дифференциальные уравнения и как их решатьСкачать