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

Метод Гира

y(x)= +y0+y01(x—x0)+y012(x—x0)(x—x1)+ y0123(x—x0)(x—x1)(x—x2)+y01234(x—x0)(x—x1)(x—x2)(x—x3),(7.25)

где y01, y012, y0123, y01234 — разделенные разности порядков с первого по четвертый.

Левую часть уравнения (7.2), т.е. производную y’(x), приближенно найдем путем дифференцирования по x полинома (7.25):

y’(x)= + —y01+y012(2x—x0—x1)+y0123[3x 2 -2x(x0+x1+x2)+ x0x1+x0x2+x1x2]+ y01234[4x 3 -3x 2 (x0+x1+x2)+2x(x0x1+x0x2+x0x3+x1x2+x1x3+x2x3)- x0x1x2—x0x1x3—x0x2x3—x1x2x3];(7.26)

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

y01=(y1-y0) / h,
y012=(y2—2y1+y0) / (2h 2 ),(7.27)
y0123=(y3-3y2+3y1— y0) / (6h 3 ),
y01234=(y4-4y3+6y2-4y1+y0) / (24h 4 ).

Полагая в (7.26) x=x4 и учитывая (7.27), получим:

y’(x4)=(3y0-16y1+36y2-48y3+25y4) / (12h).(7.28)

C другой стороны, исходное дифференциальное уравнение (7.2) при x=x4 принимает вид:

y’(x4)=f(x4,y4).(7.29)

Приравнивая правые части (7.28) и (7.29), находим:

y4=[3(4hf(x4,y4)-y0)+16y1-36y2+48y3] / 25.(7.30)

Формула (7.30) представляет собой неявную схему Гира четвертого порядка для решения задачи Коши (7.2,7.2’). Выражение (7.30) есть уравнение относительно y4, для решения которого можно применить метод простых итераций. Начальное приближение к y4 можно получить из следующих соображений. Полагая в выражении (7.26) x=x3, имеем:

y’(x3)=(-y0+6y1-18y2+10y3—y4) / (12h).(7.31)

Приравнивая правые части (7.2) при x=x3 и выражения (7.31), получим так называемую схему прогноза

y4=4hf(x3,y3)+(y0-10y3)/3-2y1+6y2,(7.32)

которую и можно использовать в качестве начального приближения для решения уравнения (7.30).

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

Рассмотрим задачу Коши для системы обыкновенных дифференциальных уравнений $$ \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

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

Пусть \( \lambda_i(t) \), \( i = 1, 2, \ldots, m \) — собственные числа матрицы $$ \begin A(t) = \< a_(t) \>, \quad a_(t) = \frac<\partial f_i><\partial u_j>(t, w). \end $$ Система уравнений (3) является жесткой, если число $$ \begin 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) область его устойчивости содержит угол $$ \begin |\arg(-\mu)| —>

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

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

Математическое описание

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

методом Гира. Решение вычисляется в одной точке XK, которая является концом интервала интегрирования.

Метод Гира для нежесткой системы является многошаговым предсказывающе — исправляющим методом Адамса, записанным в форме Нордсика, при этом предсказание и исправление имеют один и тот же порядок.

B случае, когда система уравнений является жесткой, интегрирование осуществляется специальным методом, основанном на методе типа Адамса и использующим якобиан ( ∂F/∂Y ) системы, который вычисляется подпрограммой по формулам численного дифференцирования. При интегрировании данной системы уравнений численное решение проверяется на точность; считается что значение решения в узле xn вычислено с требуемой точностью ЕРS, если выполняется следующее соотношение:

где δI — одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I — й компоненты решения в узле xn подразумевается отношение абсолютной погрешности eI этого значения в узле xn к абсолютной величине значения I — й компоненты в предыдущем узле xn — 1, т.е. eI / | yI n — 1 |, а под стандартной погрешностью — отношение eI / YPM (I), где

Tип погрешности специфицируется пользователем при обращении к подпрограмме.

Gear C.W. The automatic integration of ordinary differential equations. Communicatuons of the ACM, 14, 3 (March 1971), 176-179.

Gear C.W., Numerical Initial Value Problems in Ordinary Differential Equations, Prentice — Hall, Englewood Cliffs, N.J., 1971.

Gear C.W., The automatic integration of stiff ordinary differential equations. Information Processing 68, A.J.H.

Использование

F —имя подпрограммы вычисления значений правой части дифференциального уравнения. Первый оператоp подпрограммы должен иметь вид:
procedure F (X :Real; var Y :Array of Real; var DY :Array of Real; M :Integer);
Здесь: X, Y — значения независимой и зависимой переменных, соответственно. Вычисленное значение правой части должно быть помещено в DY. B случае системы уравнений, т.е. когда M ≠ 1, параметры Y и DY представляют массивы длины M (тип параметров X, Y и DY: вещественный);
M —количество уравнений в системе (тип: целый);
XN, YN —начальные значения аргумента и решения; в случае системы уравнений (т.е. M ≠ 1) YN представляет массив длины M (тип: вещественный);
XK —значение аргумента, при котоpом требуется вычислить решение задачи Коши (конец интервала интегрирования); XK может быть больше, меньше, или pавно XN (тип: вещественный);
HMIN —минимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений; это значение должно быть много меньше среднего ожидаемого шага интегрирования, задаваемого параметром H (тип: вещественный);
HMAX —максимальное значение абсолютной величины шага, котоpое разрешается использовать при интегрировании данной системы уравнений (тип: вещественный);
EPS —допустимая погрешность, с которой требуется вычислить все компоненты решения; тип погрешности специфицируется с помощью параметpа IU (тип: вещественный);
ISTIFJ —целый указатель метода численного интегрирования:
ISTIFJ=0 —интегрирование системы ведется методом Адамса;
ISTIFJ=1 —интегрирование ведется специальным методом, предназначенным для жестких систем;
IORDER —целая переменная, указывающая максимальный допустимый порядок метода; IORDER должен быть не больше 7 для метода Адамса и не больше 6 для метода интегрирования жестких систем;
IU —целый указатель типа погрешности численного решения:
IU = 1 —для стандартной погрешности;
IU = 2 —для относительной погрешности;
IU = 3 —для абсолютной погрешности;
H —вещественная переменная, содержащая начальное значение шага интегрирования; может задаваться с учетом направления итегрирования, т.е. положительным, если XK > XN, отрицательным, если XK < XN, или без такого учета в виде абсолютной величины;
Y —искомое решение задачи Коши, вычисленное подпрограммой для значения аргумента XK; для системы уравнений (когда M ≠ 1) задается массивом длиной M; в случае совпадения значений параметров XN и XK значение Y полагается равным начальному значению YN (тип: вещественный);
YPM —
DELTY		одномерные вещественные рабочие массивы длиной M;
RAB —одномерный вещественный рабочий массив; при интегрировании нежесткой системы уравнений RAB имеет размер 17*M, при интегрировании жесткой системы — M*(M + 17);
YP —двумерный вещественный рабочий массив размеpа 8*M;
IERR —целая переменная, служащая для сообщения об ошибках, обнаруженных в процессе работы подпрограммы; при этом:
IERR= 1 —когда неправильно задан параметр IORDER, а именно, когда IORDER превосходит максимальный допустимый порядок метода; в этом случае интегрирование системы ведется методом Гира порядка не выше 7 для нежесткой системы, и не выше 6 для жесткой;
IERR=65 —когда решение системы не может быть вычислено с требуемой точностью EPS при заданных начальном шаге H, его минимальном значении HMIN и порядке метода IORDER;
IERR=66 —когда приближенное значение решения не может быть вычислено, т.к. итерационный процесс его определения не сходится для шагов интегрирования H, больших заданного минимального значения HMIN;
IERR=67 —когда требуемая точность EPS вычисления приближенного решения меньше той, которая может быть достигнута для данной задачи при тех размерах шага интегрирования, начальное значение которого задано параметром H;
при IERR = 65, 66, 67 интегрирование системы прекращается; при желании интегрирование можно повторить обращением к подпрограмме с новыми значениями параметров H, HMIN и IORDER;
IERR=68 —когда приближенное значение решения для жесткой системы не может быть вычислено с заданной точностью; для достижения тебуемой точности следует воспользоваться версиями подпрограммы DE23E, DE25R или DE25E.
DE23E —вычисление решения задачи Коши для нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с расширенной (Extended) точностью. При этом параметры XN, YN, XK, HMIN, HMAX, EPS, H, Y, YPM, DELTY, RAB, YP и параметры X, Y и DY в подпрограмме F должны иметь тип Extended.
DE25R —вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с автоматическим выбором шага. Первый оператор подпрограммы имеет вид:
Здесь: FJ — имя подпрограммы вычисления якобиана правой части системы; первый оператор этой подпрограммы имеет вид:
procedure FJ (X :Real; var Y :Array of Real; var Z :Array of Real; M :Integer);
Здесь: X, Y — значения независимой и зависимой переменных, соответственно, причем Y представляет одномерный массив длины M; вычисленное значение якобиана должно быть помещено в двумерный массив Z размера M*M, при этом частная производная от правой части I — ого уравнения по J — ой переменной Y (J) запоминается в элементе Z (I, J) (тип параметров X, Y и Z: вещественный). Остальные параметры подпрограммы DE25R имеют тот же смысл, что и одноименные параметры подпрограммы DE23R.
DE25E —вычисление решения задачи Коши для жесткой системы обыкновенных дифференциальных уравнений первого порядка в конце интервала интегрирования методом Гира с расширенной (Extended) точностью. Первый оператор подпрограммы имеет тот же вид, что и в подпрограмме DE25R; при этом параметры XN, YN, XK, HMIN, HMAX, EPS, H, Y, YPM, DELTY, RAB, YP и параметры X, Y, DY и Z в подпрограммах F и FJ должны иметь тип Extended.
DE21R —выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы DE23R.
DE21E —выполнение одного шага численного интегрирования нежесткой и жесткой систем обыкновенных дифференциальных уравнений первого порядка методом Гира с расширенной (Extended) точностью; вызывается при работе подпрограммы DE23E.
DE24R —выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с контролем точности; вызывается при работе подпрограммы DE25R.
DE24E —выполнение одного шага численного интегрирования жесткой системы обыкновенных дифференциальных уравнений первого порядка методом Гира с расширенной (Extended) точностью; вызывается при работе подпрограммы DE25E.
UTDE12 —подпрограмма выдачи диагностических сообщений при работе подпрограмм DE21R, DE23R, DE24R, DE25R.
UTDE13 —подпрограмма выдачи диагностических сообщений при работе подпрограмм DE21E, DE23E, DE24E, DE25E.

B общем случае заданая точность EPS не гарантируется.

При работе подпрограммы и ее версий значения параметров M, XN, YN, XK, HMIN, HMAX, EPS, ISTIFJ, IORDER, IU сохраняются. Если после работы подпрограммы нет необходимости иметь начальное значение решения YN, то параметры YN и Y при обращении к ней можно совместить.


источники:

http://slemeshevsky.github.io/num-mmf/ode/html/._ode-FlatUI001.html

http://num-anal.srcc.msu.ru/lib_na/cat/de_htm_p/de23r_p.htm

© 2023 Русский язык. Привила написания. Наш первый проект тут и тут