Метод эйлера для решения дифференциальных уравнений паскаль

Метод эйлера для решения дифференциальных уравнений паскаль

Nickolay.info. Алгоритмы. Решение задачи Коши методом Эйлера

1. Типы задач для обыкновенных дифференциальных уравнений

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

В дифференциальное уравнение n-го порядка в качестве неизвестных величин входят функция у(х) и её первые n производных по аргументу х.

(1)

Из теории ОДУ известно, что уравнение (1) эквивалентно системе n уравнений первого порядка

(2)

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

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

Для системы ОДУ типа (2) начальные условия задаются в виде

(3)

Ко второму типу задач относятся так называемые граничные, или краевые задачи, в которых дополнительные условия задаются в виде функциональных соотношений между искомыми решениями. Количество условий должно совпадать с порядком n уравнения или системы. Если решение задачи определяется в интервале x принадлежит [х₀, х_k], то такие условия могут быть заданы как на границах, так и внутри интервала. Минимальный порядок ОДУ, для которых может быть сформулирована граничная задача, равен двум.

Третий тип задач для ОДУ — это задачи на собственные значения. Такие задачи отличаются тем, что кроме искомых функций у(х) и их производных в уравнения входят дополнительно m неизвестных параметров λ₁, λ₂, …, λ_m, которые называются собственными значениями. Для единственности решения на интервале [х₀, х_k] необходимо задать n+m граничных условий. В качестве примера можно назвать задачи определения собственных частот, коэффициентов диссипации, структуры электромагнитных полей и механических напряжений в колебательных системах, задачи нахождения фазовых коэффициентов, коэффициентов затухания, распределения напряженностей полей волновых процессов и т. д.

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

Большинство методов решения ОДУ основано на задаче Коши, алгоритм решения которой рассматривается далее.

Систему ОДУ (2) часто удается представить в каноническом виде, в так называемой форме Коши

(4)

При формулировке задачи Коши система (4) дополняется начальными условиями (3). Для простоты рассмотрим задачу Коши для одного уравнения типа (4), а затем полученные алгоритмы обобщим на систему n уравнений

(5)

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

(6)

который можно применить для приближенного определения искомой функции у(х). В точке х₀ + h при малых значениях h можно ограничиться двумя членами ряда (6), тогда

(7)

где O(h 2 ) — бесконечно малая величина порядка h 2 . Заменим производную у'(x₀), входящую в формулу (7), на правую часть уравнения (5):

(8)

Теперь приближенное решение в точке х₁ = х₀ + h можно вновь рассматривать как начальное условие и по формуле (8) найти значение искомой функции в следующей точке х₂ = x₁ + h₁. В результате получен простейший алгоритм решения задачи Коши, который называется методом Эйлера, или методом ломаных. Последнее название связано с геометрической интерпретацией процесса (см. рис.); искомую функцию у(х) мы заменяем ломаной линией, представляющей собой отрезки касательных к этой функции в узлах

Рис. Метод Эйлера

Формула (8) может быть получена из других соображений. Заменим производную в левой части уравнения (5) приближенным конечно-разностным отношением

Нетрудно видеть эквивалентность последнего выражения с алгоритмом Эйлера (8).

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

, (9)

где y_h(x) — приближенное решение дифференциального уравнения в точке х, полученное с шагом h; у_kh(х) — приближенное решение того же уравнения с шагом kh; р — порядок метода.

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

Для уточнения решения применима вторая формула Рунге

(10)

Формула Эйлера (8) обобщается для систем ОДУ, записанных в форме Коши (4) с начальными условиями (3)

(11)

Представленная ниже программа реализует метод Эйлера решения задачи Коши для системы дифференциальных уравнений. Функция уравнения задаётся подпрограммой f(x), точное решение – подпрограммой ft(x). Пользователь вводит интервал поиска решения [A,B], число шагов N, начальное значение Y(0). Программа выводит найденное решение и оценивает его максимальную погрешность.

Метод эйлера для решения дифференциальных уравнений паскаль

Извините — ничего не понял! Уравнение вполне себе легко решается аналитически, и выглядит это так:

1. Записываем в виде y’ — y = 1 + x.
2. Решаем однородное уравнение y’ — y = 0. Его решением, очевидно, является C*exp(x).
3. Частное решение ищем методом вариации постоянной, т.е. в виде C(x)*exp(x). Подставляя в исходное уравнение и интегрируя, находим C(x)=C1-exp(-x)*(x+2).
4. Таким образом, окончательно y=C1*exp(x)-x-2.

И при чем тут Паскаль, программирование, точность и всё такое прочее? Кроме того, полагаю, что Ваше предположение о том, что у=1 имеет какое-то отношение к таинственному «шагу», есть полный абсурд. Скорее всего, у=1 — это значение искомой функции при некотором значении аргумента х (например, при х=0), которое Вы пропустили, и которое требуется для определения константы С1.

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

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

Используя векторные обозначения, задачу (1), (2) можно записать как задачу Коши $$ \begin \tag <3>\frac> &= \pmb(t, \pmb), \quad t > 0, \\ \tag <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 \) при \( \tau\to 0 \).

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

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

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

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

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

Четвертый шаг заключается в получении численного алгоритма. Из (5) следует, что мы должны знать значение \( y^n \) для того, чтобы решить уравнение (5) относительно \( y^ \) и получить формулу для нахождения приближенного значения искомой функции на следующем временном слое \( t_ \): $$ \begin \tag <6>\pmb^ = \pmb^n + \tau \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) используется следующий метод: $$ \begin \tag <7>\frac<\pmb^ — \pmb^n> <\tau>= \pmb(t_, \pmb^). \end $$

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

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

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

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

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

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

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

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

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

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах \( \pmb^n \) и \( \pmb^ \) — один шаг по переменной \( t \). Линейный \( m \)-шаговый разностный метод записывается в виде $$ \begin \tag <13>\frac<1> <\tau>\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 \ne 0 \). Для начала расчетов по рекуррентной формуле (13) необходимо задать \( m \) начальных значений \( \pmb^0 \), \( \pmb^1 \), \( \dots \), \( \pmb^ \) (например, можно использовать для их вычисления метод Эйлера).

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

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

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

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

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

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

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

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