Явный и неявный метод эйлера

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ egin ag <1>frac &= f_i (t, u_1, u_2, ldots, u_n), quad t > 0\ ag <2>u_i(0) &= u_i^0, quad i = 1, 2, ldots, m. end $$

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ egin ag <3>frac> &= pmb(t, pmb), quad t > 0, \ ag <4>pmb(0) &= pmb_0 end $$ В задаче Коши необходимо по известному решению в точке ( t = 0 ) необходимо найти из уравнения (3) решение при других ( t ).

Существует большое количество методов численного решения задачи (3), (4). Вначале рассмотрим простейший явный метод Эйлера и его программную реализацию. Затем будут представлены методы Рунге—Кутта и многошаговые методы.

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

Идея численных методов решения задачи (3), (4) состоит из четырех частей:

1. Вводится расчетная сетка по переменной ( t ) (время) из ( N_t + 1 ) точки ( t_0 ), ( t_1 ), ( ldots ), ( t_ ). Нужно найти значения неизвестной функции ( pmb ) в узлах сетки ( t_n ). Обозначим через ( pmb^n ) приближенное значение ( pmb(t_n) ).

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

3. Аппроксимируем производные конечными разностями.

4. Формулируем алгоритм, который вычисляет новые значения ( pmb^ ) на основе предыдущих вычисленных значений ( pmb^k ), ( k 0 ) при ( au o 0 ).

Явный метод Эйлера

Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами ( t_n ) и ( t_ ): $$ omega_ au = < t_n = n au, n = 0, 1, ldots, N_t >. $$

Затем, предполагаем, что уравнение выполнено в узлах сетки, т.е.: $$ pmb^prime (t_n) = pmb(t_n, u(t_n)), quad t_n in omega_ au. $$

Заменяем производные конечными разностями. С этой целью, нам нужно знать конкретные формулы, как производные могут быть аппроксимированы конечными разностями. Простейший подход заключается в использовании определения производной: $$ pmb^prime(t) = lim_ < au o 0>frac<pmb(t+ au) — pmb(t)>< au>. $$

В произвольном узле сетки ( t_n ) это определение можно переписать в виде: $$ egin pmb^prime(t_n) = lim_ < au o 0>frac<pmb(t_n+ au) — pmb(t_n)>< au>. end $$ Вместо того, чтобы устремлять шаг сетки к нулю, мы можем использовать малый шаг ( au ), который даст численное приближение ( u^prime(t_n) ): $$ egin pmb^prime(t_n) approx frac<pmb^ — pmb^>< au>. end $$ Такая аппроксимация известна как разностная производная вперед и имеет первый порядок по ( au ), т.е. ( O( au) ). Теперь можно использовать аппроксимацию производной. Таким образом получим явный метод Эйлера: $$ egin ag <5>frac<pmb^ — pmb^n> < au>= pmb(t_n, pmb^). end $$

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение ( y^n ) для того, чтобы решить уравнение (5) относительно ( y^ ) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое ( t_ ): $$ egin ag <6>pmb^ = pmb^n + au pmb(t_n, pmb^) end $$

При условии, что у нас известно начальное значение ( pmb^0 = pmb_0 ), мы можем использовать (6) для нахождения решений на последующих временных слоях.

Программная реализация явного метода Эйлера

Выражение (6) может быть как скалярным так и векторным уравнением. И в скалярном и в векторном случае на языке Python его можно реализовать следующим образом

При решении системы (векторный случай), u[n] — одномерный массив numpy длины ( m+1 ) (( m ) — размерность задачи), а функция F должна возвращать numpy -массив размерности ( m+1 ), t[n] — значение в момент времени ( t_n ).

Таким образом численное решение на отрезке ( [0, T] ) должно быть представлено двумерным массивом, инициализируемым нулями u = np.zeros((N_t+1, m+1)) . Первый индекс соответствует временному слою, а второй компоненте вектора решения на соответствующем временном слое. Использование только одного индекса, u[n] или, что то же самое, u[n, :] , соответствует всем компонентам вектора решения.

Функция euler решения системы уравнений реализована в файле euler.py:

Строка F_ = lambda . требует пояснений. Для пользователя, решающего систему ОДУ, удобно задавать функцию правой части в виде списка компонент. Можно, конечно, требовать чтобы пользователь возвращал из функции массив numpy , но очень легко осуществлять преобразование в самой функции решателе. Чтобы быть уверенным, что результат F будет нужным массивом, который можно использовать в векторных вычислениях, мы вводим новую функцию F_ , которая вызывает пользовательскую функцию F «прогоняет» результат через функцию assaray модуля numpy .

Неявный метод Эйлера

При построении неявного метода Эйлера значение функции ( F ) берется на новом временном слое, т.е. для решении задачи (5) используется следующий метод: $$ egin ag <7>frac<pmb^ — pmb^n> < au>= pmb(t_, pmb^). end $$

Таким образом для нахождения приближенного значения искомой функции на новом временном слое ( t_ ) нужно решить нелинейное уравнение относительно ( pmb^ ): $$ egin ag <8>pmb^ — au pmb(t_, pmb^) — y^n = 0. end $$

Для решения уравнения (8) можно использовать, например, метод Ньютона.

Программная реализация неявного метода Эйлера

Функция backward_euler решения системы уравнений реализована в файле euler.py:

Отметим, что для нахождения значения u[n+1] используется функция fsolve модуля optimize библиотеки scipy . В качестве начального приближения для решения нелинейного уравнения используется значение искомой функции с предыдущего слоя u[n] .

Методы Рунге—Кутта

Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ egin ag <9>frac<pmb^ — pmb^n> < au>= sum_^s b_i pmb_i, end $$ где $$ egin ag <10>pmb_i = pmbleft( t_n + c_i au, pmb^n + au sum_^s a_pmb_j
ight), quad i = 1, 2, ldots, s. end $$ Формула (9) основана на ( s ) вычислениях функции ( pmb ) и называется ( s )-стадийной. Если ( a_ = 0 ) при ( j geq i ) имеем явный метод Рунге—Кутта. Если ( a_ = 0 ) при ( j > i ) и ( a_
e 0 ), то ( pmb_i ) определяется неявно из уравнения $$ egin ag <11>pmb_i = pmbleft( t_n + c_i au, pmb^n + au sum_^ a_pmb_j + au a_ pmb_i
ight), quad i = 1, 2, ldots, s. end $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.

Одним из наиболее распространенных является явный метод Рунге-Кутта четвертого порядка: $$ egin ag <12>pmb_1 & = pmb(t_n, pmb^n), &quad pmb_2 &= pmbleft( t_n + frac< au><2>, pmb^n + au frac<pmb_1> <2>
ight),\ pmb_3 &= pmbleft( t_n + frac< au><2>, pmb^n + au frac<pmb_2> <2>
ight), &quad pmb_4 &= pmbleft( t_n + au, pmb^n + au pmb_3
ight),\ frac<pmb^ -pmb^n> < au>&= frac<1> <6>(pmb_1 + 2pmb_2 + 2pmb_3 + pmb_4) & & end $$

Многошаговые методы

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах ( pmb^n ) и ( pmb^ ) — один шаг по переменной ( t ). Линейный ( m )-шаговый разностный метод записывается в виде $$ egin ag <13>frac<1> < au>sum_^m a_i pmb^ = sum_^ b_i pmb(t_, pmb^), quad n = m-1, m, ldots end $$ Вариант численного метода определяется заданием коэффициентов ( a_i ), ( b_i ), ( i = 0, 1, ldots, m ), причем ( a_0
e 0 ). Для начала расчетов по рекуррентной формуле (13) необходимо задать ( m ) начальных значений ( pmb^0 ), ( pmb^1 ), ( dots ), ( pmb^ ) (например, можно использовать для их вычисления метод Эйлера).

Различные варианты многошаговых методов (методы Адамса) решения задачи с начальными условиями для систем обыкновенных дифференциальных уравнений могут быть получены на основе использования квадратурных формул для правой части равенства $$ egin ag <14>pmb(t_) — pmb(t_n) = int_^> pmb(t, pmb) dt end $$

Для получения неявного многошагового метода используем для подынтегральной функции интерполяционную формулу по значениям функции ( pmb^ = pmb(t_, pmb^) ), ( pmb^n ), ( dots ), ( pmb^ ), т.е. $$ egin ag <15>frac<pmb^ — pmb^n> < au>= sum_^ b_i pmb(t_, pmb^) end $$

Для интерполяционного метода Адамса (15) наивысший порядок аппроксимации равен ( m+1 ).

Для построения явных многошаговых методов можно использовать процедуру экстраполяции подынтегральной функции в правой части (14). В этом случае приближение осуществляется по значениям ( pmb^n ), ( pmb^ ), ( dots ), ( pmb^ ) и поэтому $$ egin ag <16>frac<pmb^ — pmb^n> < au>= sum_^ b_i pmb(t_, pmb^) end $$

Для экстраполяционного метода Адамса (16) погрешность аппроксимации имеет ( m )-ый порядок.

На основе методов Адамса строятся и схемы предиктор–корректор. На этапе предиктор используется явный метод Адамса, на этапе корректора — аналог неявного метода Адамса. Например, при использовании методов третьего порядка аппроксимации в соответствии с (18) для предсказания решения положим $$ frac<pmb^ — pmb^n> < au>= frac<1> <12>(23 pmb^ -16pmb^ + 5pmb^). $$ Для уточнеия решения (см. (17)) используется схема $$ frac<pmb^ — pmb^n> < au>= frac<1> <24>(9pmb^ + 19pmb^ — 5pmb^ + pmb^). $$ Аналогично строятся и другие классы многошаговых методов.

Жесткие системы ОДУ

При численном решении задачи Коши для систем обыкновенных дифференциальных уравнений (3), (4) могут возникнуть дополнительные трудности, порожденные жесткостью системы. Локальные особенности поведения решения в точке ( u = w ) передаются линейной системой $$ egin frac

= sum_^ frac<partial f_i> <partial u_j>(t, w) v + ar(t), quad t > 0. end $$

Пусть ( lambda_i(t) ), ( i = 1, 2, ldots, m ) — собственные числа матрицы $$ egin A(t) = < a_(t) >, quad a_(t) = frac<partial f_i><partial u_j>(t, w). end $$ Система уравнений (3) является жесткой, если число $$ egin S(t) = frac <max_<1 leq i leq m>|Re lambda_i(t)|> <min_<1 leq i leq m>|Re lambda_i(t)|> end $$ велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной ( t ).

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

Метод называется ( A )-устойчивым, если при решении задачи Коши для системы (3) область его устойчивости содержит угол $$ egin |arg(-mu)| —>

Дата добавления: 2015-06-12 ; просмотров: 10685 ; Нарушение авторских прав

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

Если задано дифференциальное уравнение

(3.1)

и начальные условия , то очередное значение может быть получено интегрированием (3.1):

(3.2)

Определенный интеграл в (3.2) численно равен площади под кривой на интервале (рис. 3.2).

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

Подставив в (3.2) приближенные значения интеграла, можно получить две формулы:

(3.3)

. (3.4)

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

Формула (3.4) соответствует неявному методу Эйлера. Здесь в правой части выражения используется неизвестное значение , поэтому вычислить его непосредственно по этой формуле нельзя.

Более точное значение интеграла (3.2) дает метод трапеций, которому соответствует отрезок 3 на рис. 3.3. Тогда

. (3.5)

Эта формула относится, очевидно, тоже к неявным.

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

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

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

Пусть задано нелинейное алгебраическое уравнение

(3.6)

где – вектор переменных.

Разложение (3.6) в ряд Тэйлора с сохранением только линейных членов дает приближенную замену

(3.7)

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

– неизвестное значение переменной на шаге интегрирования.

Выражение (3.7) может быть записано как линейное алгебраическое уравнение

,(3.8)

где – вычисляется для известных значений переменных на предыдущем шаге интегрирования;

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

, (3.9)

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

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

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

Явные формулы метода Эйлера имеют вид

Неявные формулы запишутся следующим образом:

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

Здесь ,

Матричная запись имеет вид

.

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

3.2.3. Выбор между явными и неявными методами
в процедурах моделирования технических систем

Выбор между явными и неявными методами представляет серьезную проблему. Многие специалисты считают неявные методы более мощным и универсальным инструментом для решения задач моделирования технических систем [23, 15]. Следует, однако, заметить, что лишь недавно появились достаточно мощные и универсальные системы автоматизированного моделирования, такие, как, например, MATLAB или МВТУ [17], допускающие выбор явного или неявного метода решения задачи. Раньше использовались либо явные, либо неявные методы, так как это требовало разных компонентных моделей.

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

Рассмотрим эту особенность неявных методов на примере явного и неявного методов Эйлера [21], определяемых формулами (3.3) и (3.4), соответственно.

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

.

Характеристическое уравнение данной динамической системы имеет вид

, или ,

где – постоянная времени системы.

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

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

.

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

или .

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

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

Переход к методам более высокого порядка мало меняет картину. Для метода Рунге – Кутты 4-го порядка требование устойчивости ограничивает шаг величиной , или, в более общем виде, , где – максимальное собственное значение матрицы Якоби [29].

Применение неявного метода Эйлера к той же системе дает

,

где ограничение на величину шага выглядит по-другому:

,

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

Формой Коши для системы ОДУ является форма, в которой уравнения системы разрешены относительно производных:

(4.7)

Система (4.7) для решения задачи Коши дополняется начальными условиями. В пособии для упрощения выкладок будем решать задачу Коши для одного уравнения:

. (4.8)

Из теории рядов известно, что любую функцию y(x) в малой окрестности точки х₀ можно разложить в ряд по степеням х, который носит название ряда Тейлора:

(4.9)

Ограничив ряд при малом шаге двумя членами, получим

(4.10)

где O(h 2 ) — бесконечно малая величина порядка h 2 .

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

(4.11)

Рассчитав значение функции в точке , принимаем его за начальное для следующего шага и повторяем расчет. Для любого i-го шага интегрирования будем иметь:

или (4.12)

Такая формула носит название рекуррентной. Метод численного интегрирования, использующий два члена разложения ряда Тейлора, носит название явного метода Эйлера. Явным этот и другие подобные методы называются потому, что в уравнении используется уже известное значение производной — в точке x_i , то есть производ

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

Как известно, уравнение прямой имеет вид y=ax+b, где а — тангенс угла между осью абцисс и прямой, b — значение y при x, равном нулю. Также известно, что геометрический смысл производной — величина производной f’(x) равна тангенсу угла наклона касательной, проведенной к графику кривой функции в данной точке.

Принимая за аргумент шаг h , можно построить на графике прямую линию из точки x₀,y₀ в точку x₁,y₁ под углом, тангенс которого равен производной в точке x₀ — f(x₀,y₀) (рис.4.2).

Из точки х₁ под углом, тангенс которого равен f(x₁,y₁₎, проводим прямую в точку х₂ и т. д. В качестве решения получим ломаную линию. В связи с этим метод Эйлера еще носит название метода ломаных.

Подставив в формулу (4.12) вместо производной в исходной точке f(x_i,y_i) значение еще неизвестной производной в искомой точке x_i+h — f(x_i+1,y_i+1), получим формулу неявного метода Эйлера:

(4.13)

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

, (4.14)

что в задачах электротехники, как правило, имеет место (см. уравнение (1.3) в главе 1 — ), от итерационного метода можно перейти к прямому методу интегрирования неявным методом Эйлера. Подставив выражение для производной из (4.14) в (4.13), получим:

.

Разрешив данное уравнение относительно y_{i+1 ,}получим прямую формулу для интегрирования ОДУ неявным методом Эйлера:

. (4.15)

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

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

На каждом шаге метода Эйлера решение y(x) определяется с погрешностью за счет отбрасывания членов ряда Тейлора, пропорциональных h в степени выше первой. Это означает, что локальная погрешность (погрешность на шаге интегрирования) метода Эйлера
имеет второй порядок. Глобальная погрешность (погрешность на всем интервале интегрирования) имеет первый порядок.

При постоянном шаге h для оценки погрешности применима формула Рунге

, (4.16)

где y_h(x) — приближенное решение ОДУ в точке х, полученное с шагом h; y_kh(x) — приближенное решение того же уравнения, полученное с шагом kh; p — порядок метода.

Процедуры, реализующие явный и неявный методы Эйлера, приведены в конце главы (ПРОГРАММЫ 4.1 и 4.2). Сами процедуры дают решение на шаге интегрирования. Для поиска решения на всем интервале изменения аргумента в вызывающей программе должен быть организовано циклическое изменение аргумента х.

Как видно, метод Эйлера как в явной, так и в неявной форме достаточно просто реализуется на ЭВМ. Коэффициенты А и В неявного метода рассчитываются в вызывающей программе.

Для повышения точности решения разработаны одношаговые методы, позволяющие учесть большее количество членов ряда Тейлора, чем в методе Эйлера. Методов достаточно много [9]. Рассмотрим некоторые из них.

Не нашли то, что искали? Воспользуйтесь поиском: