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

где y₀₁, y₀₁₂, y₀₁₂₃, y₀₁₂₃₄ — разделенные разности порядков с первого по четвертый.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Проиллюстрируем указанные шаги. Для начала введем расчетную сетку. Очень часто сетка является равномерной, т.е. имеет одинаковое расстояние между узлами ( t_n ) и ( t_ ): $$ omega_tau = . $$

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Одношаговый метод Рунге—Кутта в общем виде записывается следующим образом: $$ begin tag frac<pmb^ — pmb^n> = sum_^s b_i pmb_i, end $$ где $$ begin tag pmb_i = pmbleft( t_n + c_itau, 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 pmb_i = pmbleft( t_n + c_itau, pmb^n + tau sum_^ a_pmb_j + tau a_ pmb_i right), quad i = 1, 2, ldots, s. end $$ О таком методе Рунге—Кутта говорят как о диагонально-неявном.

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

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

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

В методах Рунге—Кутта в вычислениях участвуют значения приближенного решения только в двух соседних узлах ( pmb^n ) и ( pmb^ ) — один шаг по переменной ( t ). Линейный ( m )-шаговый разностный метод записывается в виде $$ begin tag frac 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 pmb(t_) — pmb(t_n) = int_^<t_> pmb(t, pmb) dt end $$

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

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

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

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

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

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

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

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

= sum_^ frac (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(t, w). end $$ Система уравнений (3) является жесткой, если число $$ begin S(t) = frac <max_|Re lambda_i(t)|> <min_|Re lambda_i(t)|> end $$ велико. Это означает, что в решении присутствуют составляющие с сильно различающимися масштабами изменения по переменной ( t ).

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

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

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

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

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

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

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

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

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

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

где δ_I — одна из погрешностей следующих типов: абсолютная, относительная или стандартная. При этом под относительной погрешностью приближенного значения I — й компоненты решения в узле x_n подразумевается отношение абсолютной погрешности e_I этого значения в узле x_n к абсолютной величине значения I — й компоненты в предыдущем узле x_{n — 1}, т.е. e_I / | y_I n — 1 |, а под стандартной погрешностью — отношение e_I / 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 при обращении к ней можно совместить. 🎬 Видео Операторный метод решения дифференциальных уравнений | Решение задачСкачать Численное решение дифференциальных уравнений (задачи Коши)Скачать Видеоурок "Системы дифференциальных уравнений"Скачать Численное решение задачи Коши методом ЭйлераСкачать Линейное дифференциальное уравнение Коши-ЭйлераСкачать Методы Гира для решения жестких СОДУ. Gear's Methods for Stiff Systems of ODE.Скачать Видеоурок "Нахождение частных решений по виду правой части"Скачать Поделиться или сохранить к себе: Смотрите также Сера — химические свойства, получение, соединения. Нитрат кальция: способы получения и химические свойства Кальций: способы получения и химические свойства Гидроксид натрия: способы получения и химические свойства Гидроксид кальция: способы получения и химические свойства Получение бензола Разница между глюкозой и сахарозой Гидроксид калия: способы получения и химические свойства