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

Программа DE10 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом Рунге-Кутта 4-го порядка аппроксимации с постоянным шагом без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 5 порядка включительно.
Решение системы вычисляется по формулам:

call DE10(DS, A, B, N, Y)

DS, A, B, N — входные параметры;
Y — входной и выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных в точке X :

Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . N должно быть >=1;
Real Y(1:m) — массив, который на входе должен содержать начальные значения функций в точке A , на выходе будет содержать решение системы в точке B .

Программа DE15 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом Рунге-Кутта 5-го порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:

На каждом шаге интегрирования системы программа вычисляет асимптотическую оценку погрешности решения:

Если погрешность оказывается меньше заданного значения ε , программа переходит к следующему шагу, если больше — программа повторяет вычисления с уменьшенной величиной шага. Величина нового шага h_new расчитывается по формуле:
_{5 ______}
h_new = h_old•√ ε/e .

call DE15(DS, X, Xout, H, Eps, Y, Error)

DS, Xout — входные параметры;
X, H, Eps, Y — входные и выходные параметры;
Error — выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных (см. прог. DE10 );
Real X, Xout — начальная и конечная точки интегрирования;
Real H — начальный шаг интегрирования, с которым программе предлагается начать вычисления. Допустимо задать абсолютное значение, даже если шаг должен быть отрицательным. На выходе переменная содержит величину шага в конце промежутка интегрирования;
Real Eps — требуемая точность интегрирования на каждом шаге;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке X , на выходе будет содержать решение системы в точке Xout . Число M задает размерность решаемой системы уравнений;
Integer Error — индикатор ошибки.
Error=0 , если интегрирование системы успешно закончено и решение вычислено с заданной точностью Eps на каждом шаге;
Error=1 , если интегрирование нельзя начать, т.к. значение |Xout-X| меньше минимально-допустимой величины шага h_min , которая вычисляется программой. Параметр H приобретает значение h_min ;
Error=2 , если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины ε_min , которая вычисляется программой. Параметр Eps приобретает значение ε_min ;
Error=65 , если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps . Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.

Программа DE19 вычисляет значение производной f₀ и разности назад ∇f₀ , ∇ 2 f₀ , . ∇ k f₀ порядка k в точке x₀ (т.н. фронт метода), необходимые для начала интегрирования системы дифференциальных уравнений первого порядка по методу Адамса. Для нахождения этих разностей применяется алгоритм разгона, который основан на итерационном применении явных формул Адамса с последовательно увеличивающимся порядком аппроксимации. Используя начальное значение y₀ в точке x₀ и формулу Адамса первого порядка, вычисляем

Находим ∇f₁ (1) = f₁ (1) — f₀ и полагаем ∇f₀ = ∇f₁ (1) . Используя формулу Адамса второго порядка, вычисляем

Находим f₂ (2) , ∇f₂ (2) , ∇ 2 f₂ (2) и полагаем

Используя формулу Адамса третьего порядка, вычисляем

Последовательно применяя формулы Адамса все более высокого порядка получаем все нужные нам разности ∇f₀ , ∇ 2 f₀ , . ∇ k f₀ до заданного порядка k включительно. Более подробную информацию об алгоритме разгона можно найти в книге [Б6].

call DE19(DS, X0, Y0, H, M, Order, DF, XH, YH, Err)

DS, X0, Y0, H, M, Order — входные параметры;
DF, XH, YH, Err — выходные параметры;

DS(X, Y, DY) — процедура, которая вычисляет значения производных (см. программу DE10 );
Real X0 — начальная точка интегрирования;
Real Y0(1:M) — начальные значения функций в точке X0 ;
Real H — шаг интегрирования;
Integer M — размерность системы дифференциальных уравнений;
Integer Order — порядок метода, который будет применен для решения системы. Значение Order должно лежть в пределах от 1 до 6;
Real DF(1:M,0:Order) — массив содержит вычисленные программой значения разностей назад в точке X0 до порядка Order включительно;
Real XH — значение аргумента, равное Order*H ;
Real YH(1:M) — вычисленное програмой значение функции в точке XH ;
Real Err — асимптотическая оценка погрешности вычисленного значаения YH .

Программа DE20 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом Адамса с порядком аппроксимации от 1 до 6 без контроля точности. Предполагается, что решение системы имеет непрерывные производные на один порядок больше применяемого порядка аппроксимации.
В каждом узле интегрирования решение отыскивается по схеме прогноз-коррекция, которая заключается в следующем. Исходя из точки x_n , в которой задано значение функции y_n и фронт метода f_n , ∇f_n , ∇ 2 f_n , . ∇ k-1 f_n , вычисляется решение в точке x_n+1=x_n+h по явной формуле Адамса (прогноз):
_k-1
y_n+1 (p) = y_n + h•∑ α_j∇ j f_n , (1)
j=0
где ∇ j обозначает разность назад j -го порядка, k — порядок аппроксимации. По предсказанному значению y_n+1 (p) в точке x_n+1 вычисляется значение производной f_n+1 = f(x_n+1, y_n+1 (p) ) и новый фронт метода ∇f_n+1 , ∇ 2 f_n+1 , . ∇ k-1 f_n+1 . Теперь предсказанное значение y_n+1 (p) можно уточнить по неявной формуле Адамса (коррекция):
_k-1
y_n+1 (c) = y_n + h•∑ β_j∇ j f_n+1 . (2)
j=0
Полученное скорректированное значение принимается за значение функции в данном узле и программа переходит к вычислениям в следующем узле и так продолжается, пока не будет достигнута конечная точка интегрирования.
Итерационный процесс (1),(2) может быть преобразован к виду [Б6]:
_k-1
y_n+1 (c) = y_n+1 (p) — α_k-1•h•∑ ∇ j f_n + α_k-1•h•f_n+1 .
j=0
Коэффициенты α_k и β_k для явной и неявной формул Адамса в разностной форме для разного порядка аппроксимации k приведены в таблице 1 ниже. Увеличение поядка аппроксимации осуществляется простым добавлением новых слогаемых в формулах (1) и (2). Локальная погрешность равна коэффициенту при первом отброшенном члене суммы. В таблице 2 приведены формулы Адамса в классической форме.

call DE20(DS, A, B, N, M, Order, Y)

DS, A, B, N, M, Order — входные параметры;
Y — входной и выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных (см. программу DE10 );
Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . Значение N должно быть >= Order ;
Integer M — размерность системы дифференциальных уравнений;
Integer Order — порядок метода, который следует применить для решения системы. Значение Order должно лежть в пределах от 1 до 6;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке A , на выходе будет содержать решение системы в точке B .
В процессе вычислений программа DE20 вызывает программу DE19 .

Программа DE21 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом Адамса четвертого порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до пятого порядка включительно.
Решение системы осуществляется по схеме прогноз-коррекция, которая изложена ранее в описании к программе DE20 .
В каждом узле интегрирования программа вычисляет погрешность полученного решения. Если погрешность оказалась меньше заданного значения ε , программа переходит к следующему узлу. Если погрешность больше ε , программа уменьшает шаг в два раза и повторяет вычисления в данном узле с уменьшенной величиной шага. В том случае, когда погрешность решения оказывается меньше, чем ε/32 несколько шагов подряд, шаг интегрирования удваивается. Вычисления продолжаются, пока не будет достигнута конечная точка промежутка.
Погрешность вычисляется на основе анализа асимптотических формул остаточных членов явного и неявного методов Адамса.
ρ (p) = y (e) — y (p) = (251/720)·h 5 ·y (v) (ξ) ≈ (251/720)·h·∇ 4 f_n+1
ρ (c) = y (e) — y (c) = -(19/720)·h 5 ·y (v) (ξ) ≈ -(19/720)·h·∇ 4 f_n+1
где ρ (p) и ρ (c) — погрешности явного и неявного методов, y (e) — точное решение (неизвестное), ξ — точка из промежутка ( x_n , x_n+1 ). Вычитая из первой формулы вторую, имеем:
y (c) — y (p) = (251/720 + 19/720)·h 5 ·y (v) (ξ) = (270/720)·h 5 ·y (v) (ξ) ,
следовательно h 5 ·y (v) (ξ) = (720/270)·(y (c) — y (p) ) . Тогда
ρ (c) = -(19/720)·(720/270)·(y (c) — y (p) ) = -(19/270)·(y (c) — y (p) ) =
= -(19/270)·(3/8)·h·∇ 4 f_n+1 = -(19/720)·h·∇ 4 f_n+1
Таким образом величина погрешности совпадает с остаточным членом неявной формулы Адамса.

call DE21(DS, X, Xout, H, Eps, Y, Error)

DS, Xout — входные параметры;
X, H, Eps, Y — входные и выходные параметры;
Error — выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных (см. программу DE10 );
Real X, Xout — начальная и конечная точки интегрирования;
Real H — начальный шаг интегрирования, с которым программе предлагается начать вычисления. Допустимо задать абсолютное значение, даже если шаг должен быть отрицательным. На выходе переменная содержит величину шага в конце промежутка интегрирования;
Real Eps — требуемая точность интегрирования на каждом шаге;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке X , на выходе будет содержать решение системы в точке Xout . Число M задает размерность решаемой системы уравнений;
Integer Error — индикатор ошибки.
Error=0 , если интегрирование системы успешно закончено и решение вычислено с заданной точностью Eps на каждом шаге;
Error=1 , если интегрирование нельзя начать, т.к. значение |Xout-X| меньше минимально-допустимой величины шага h_min , которая вычисляется программой. Параметр H приобретает значение h_min ;
Error=2 , если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины ε_min , которая вычисляется программой. Параметр Eps приобретает значение ε_min ;
Error=3 , если интегрирование нельзя начать, т.к. погрешность решения, определенная программой при вычислении начальных значений, меньше Eps . Чтобы начать интегрирование, нужно или увеличить величину шага H , или уменьшить значение Eps ;
Error=65 , если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps . Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
В процессе вычислений программа DE21 вызывает программу DE19 .

Программа DE30 решает систему из m обыкновенных дифференциальных уравнений второго порядка y_i«(x) = f_i(x, y₁(x), . y_m(x), y₁‘(x), . y_m‘(x)) , i=1, 2, . m методом Рунге-Кутта 4-го порядка аппроксимации без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:

call DE30(DS2, A, B, N, Y, DY)

DS2, A, B, N — входные параметры;
Y, DY — входной и выходной параметр;

DS2(X, Y, DY, DDY) — процедура, которая вычисляет значения вторых производных в точке X :

Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . N должно быть >=1;
Real Y(1:m) — массив, который на входе должен содержать начальные значения функций в точке A , на выходе будет содержать решение системы в точке B ;
Real DY(1:m) — массив, который на входе должен содержать начальные значения производных в точке A , на выходе будет содержать значения производных в точке B .

Программа DE34 вычисляет значение производной f₀ , вспомогательной переменной z₀ и разности назад ∇f₀ , ∇ 2 f₀ , . ∇ k f₀ порядка k в точке x₀ , необходимые для начала интегрирования системы дифференциальных уравнений второго порядка по методу Штермера. Указанные величины нходятся с помощью итерационного процесса, в котором применяются явная формула Адамса (для нахождения значения производной y₀‘ ) и аналог фрмулы Адамса для уравнения второго порядка (для нахождения значения переменной z₀ ):
_k-1
y₁‘ = y₀‘ + h•∑ α_j∇ j f_n ,
j=0
_k-1
z₁ = y₀ + hy₀‘ + h•∑ μ_j∇ j f_n ,
j=0
y₁ = y₀ + hz₁ .
где α_j и μ_j — коэффициенты этих формул, j=0, 1, . k . Применяя эти формулы раз за разом c последовательно увеличивающимся порядком к исходным данным x₀ , y₀ , y₀ ‘, мы получаем необходимые значения для начала интегрирования системы по методу Штермера из точки x₀ . Более подробную информацию об этой процедуре можно найти в книге [Б6].

call DE34(DS2, X0, Y0, DY0, H, M, Order, Z, DF, XH, YH, DYH, Err)

DS2, X0, Y0, DY0, H, M, Order — входные параметры;
Z, DF, XH, YH, DYH, Err — выходные параметры;

DS2(X, Y, DY, DDY) — процедура, которая вычисляет значения вторых производных (см. программу DE30 );
Real X0 — начальная точка интегрирования;
Real Y0(1:M), DY0(1:M) — начальные значения функции и ее первой производной в точке X0 ;
Real H — шаг интегрирования;
Integer M — размерность системы дифференциальных уравнений;
Integer Order — порядок метода, который будет применен для решения системы. Значение Order должно лежть в пределах от 1 до 6;
Real Z — вычисленное программой значение вспомогательной переменной в точке X0 ;
Real DF(1:M,0:Order) — массив содержит вычисленные программой значения разностей назад в точке X0 до порядка Order включительно;
Real XH — значение аргумента, равное Order*H ;
Real YH(1:M), DYH(1:M) — вычисленное програмой значение функции и ее первой производной в точке XH ;
Real Err — асимптотическая оценка погрешности вычисленного значаения YH .

Программа DE35 решает систему из m обыкновенных дифференциальных уравнений второго порядка y_i«(x) = f_i(x, y₁(x), . y_m(x), y₁‘(x), . y_m‘(x)) , i=1, 2, . m методом Штермера с порядком аппроксимации от 1 до 6 без контроля точности. Предполагается, что решение системы имеет непрерывные производные на два порядка больше применяемого порядка аппроксимации.
В каждом узле интегрирования решение отыскивается по схеме прогноз-коррекция, которая заключается в следующем. В точке x_n+1=x_n+h вычисляется приближенное решение по явной формуле Штермера (прогноз):
_k-1
y_n+1 (p) = y_n + ∇y_n + h 2 •∑ γ_j∇ j f_n ,
j=0
где ∇ j обозначает разность назад j -го порядка, k — порядок аппроксимации. Полученное приближение можно уточнить, выполнив вычисления по неявной формуле Штермера (коррекция):
_k-1
y_n+1 (с) = y_n + ∇y_n + h 2 •∑ δ_j∇ j f_n+1 .
j=0
Здесь γ_j и δ_j обозначают коэффициенты явной и неявной формул Штермера. С целью уменьшения вычислительной погрешности при интегрировании системы вводится новая переменная z_n и формулы Штермера пиобретают вид:
_k-1
z_n+1 (p) = z_n + h•∑ γ_j∇ j f_n ,
j=0
_k-1
z_n+1 (c) = z_n + h•∑ δ_j∇ j f_n+1 ,
j=0
y_n+1 (c) = y_n + hz_n+1 (c) .
Значения первой производной, также как и в программе DE20 , вычисляются по схеме прогноз-коррекция по явной и неявной формулам Адамса:
_k-1
y_n+1‘ (p) = y_n‘ + h•∑ α_j∇ j f_n ,
j=0
_k-1
y_n+1‘ (c) = y_n‘ + h•∑ β_j∇ j f_n+1 .
j=0
Для упрощения вычислений формулы коррекции слегка видоизменяются, чтобы в их состав входили уже вычисленные предсказанные значения z_n+1 (p) и y_n+1‘ (p) . Окончательно вычислительный процесс выглядит следующим образом. Сначала по явным формулам вычисляются значения функции и производной в точке x_n+1 :
_k-1
z_n+1 (p) = z_n + h•∑ γ_j∇ j f_n ,
j=0
y_n+1 (p) = y_n + hz_n+1 (p) ,
_k-1
y_n+1‘ (p) = y_n‘ + h•∑ α_j∇ j f_n .
j=0
Далее вычисляется значение f_n+1=f(x_n+1, y_n+1 (p) , y_n+1‘ (p) ) и новый фронт метода. Скорректированные значения вычисляются по формулам:
_k-1
z_n+1 (c) = z_n+1 (p) — γ_k-1•h•∑ ∇ j f_n + γ_k-1•h•f_n+1 .
j=0
y_n+1 (c) = y_n + hz_n+1 (c) ,
_k-1
y_n+1‘ (c) = y_n+1‘ (p) — α_k-1•h•∑ ∇ j f_n + α_k-1•h•f_n+1 .
j=0
Полученные скорректированные значения принимаются за значение функции и ее производной в данном узле и программа переходит к вычислениям в следующем узле и так продолжается, пока не будет достигнута конечная точка интегрирования.
Коэффициенты γ_k и δ_k для явной и неявной формул Штермера в разностной форме для разного порядка аппроксимации k можно найти в книге [Б6].

call DE35(DS2, A, B, N, M, Order, Y, DY)

DS2, A, B, N, M, Order — входные параметры;
Y, DY — выходные параметры;

DS2(X, Y, DY, DDY) — процедура, которая вычисляет значения вторых производных (см. программу DE30 );
Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . Значение N должно быть >= Order ;
Integer M — размерность системы дифференциальных уравнений;
Integer Order — порядок метода, который следует применить для решения системы. Значение Order должно лежть в пределах от 1 до 6;
Real Y(1:M), DY(1:M) — массивы, которые на входе должны содержать начальные значения функций и их первые производные в точке A , на выходе будет содержать решение системы и значения производных в точке B .
В процессе вычислений программа DE35 вызывает программу DE34 .

Программа DE23 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом трапеций (неявным методом Адамса 1 порядка аппроксимации):
y_n+1‘ = y_n + (1/2)h(f(x_n+1, y_n+1) + f(x_n, y_n))
Для повышения устойчивости итерационный процесс для нахождения y_n+1(x) выполняется по методу Ньютона.
На каждом шаге выполняется от 1 до 3 итераций. Итерации прерываются, если каждая из компонент вектоа решения становится меньше заданного числа Eps . Если после выполнения 3-х итераций какая-либо из компонент не достигла заданной точности, специальная переменная nErr увеличивается на 1. На выходе из программы эта переменная содержит общее число шагов, на которых не была достигнута заданная точность Eps .
Для реализации итераций по методу Ньютона на каждой итерации требуется обращение матрицы Якоби, что для сложных систем уравнений может быть достаточно трудоёмко. Если от итерации к итерации матрица Якоби меняется незначительно, заданием специального значения в переменной nWex можно оставить обращение матрицы только перед первой итерацией, а в остальных итерациях пользоваться уже вычисленным значением. Это позволяет значительно уменьшить объём вычислений. В некоторых случаях (для систем уравнений с постоянными коэффициентами) матрица Якоби не зависит от аргумента и функций, т.е. является числовой матрицей. Тогда матрица Якоби обращается только один раз в начале выполнения программы и далее на всём выполнении программы используется уже вычисленное значение.
В процессе вычислений программа DE23 вызывает программу AIG2R обращения матрицы.

call DE23(DS, DSJ, A, B, Y, N, Eps, met, nWex, nErr)

DS, DSJ, A, B, N, met, nWex — входные параметры;
Y, Eps — входные и выходные параметры;
nErr — выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных;
DSJ(X, Y, DYJ) — процедура, которая вычисляет значения матрицы Якоби. Частная производня от правой части i-го уравнения по j-ой переменной Y(i) запоминается в элементе DYJ(i,j) ;
Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . N должно быть >=1;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке A , на выходе будет содержать решение системы в точке B . Число M задает размерность решаемой системы уравнений;
Real Eps — требуемая точность интегрирования на каждом шаге; если задано значение меньшее, чем минимально-допустимое, программа увеличит его до минимально-допусимого;
Character met — метод решения уавнения: met = ‘T’ — метод трапеций, met = ‘E’ — неявный метод Эйлера;
Integer nWex — задаёт режим вычисления Якобиана, может принимать значения 0, 1 или 2; 0 — Якобиан вычисляется только один раз в начале программы; 1 — вычисляется 1 раз на кждом шаге перед первой итерацией; 2 — вычисляется перед каждой итерацией, т.е. от 1 до 3 раз на каждом шаге.
Integer nErr — счётчик количества точек, при которых за 3 итерации не удалось достигнуть заданной точности по какой-либо компоненте решения.

Программа DE23 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(x, y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m по формуле дифференцирования назад 2-го поядка аппроксимации:
y_n+2‘ = (4/3)y_n+1 — (1/3)y_n + (2/3)hf(x_n+2, y_n+2)
Для повышения устойчивости итерационный процесс для нахождения y_n+2(x) выполняется по методу Ньютона. Первый шаг программа выполняет неявным методом Эйлера. На следующих шагах предсказание начального приближения осуществляется вычислением экстраполяционного многочлена Эрмита.
На каждом шаге выполняется от 1 до 3 итераций. Итерации прерываются, если каждая из компонент вектоа решения становится меньше заданного числа Eps . Если после выполнения 3-х итераций какая-либо из компонент не достигла заданной точности, специальная переменная nErr увеличивается на 1. На выходе из программы эта переменная содержит общее число шагов, на которых не была достигнута заданная точность Eps .
Для реализации итераций по методу Ньютона на каждой итерации требуется обращение матрицы Якоби, что для сложных систем уравнений может быть достаточно трудоёмко. Если от итерации к итерации матрица Якоби меняется незначительно, заданием специального значения в переменной nWex можно оставить обращение матрицы только перед первой итерацией, а в остальных итерациях пользоваться уже вычисленным значением. Это позволяет значительно уменьшить объём вычислений. В некоторых случаях (для систем уравнений с постоянными коэффициентами) матрица Якоби не зависит от аргумента и функций, т.е. является числовой матрицей. Тогда матрица Якоби обращается только один раз в начале выполнения программы и далее на всём выполнении программы используется уже вычисленное значение.
В процессе вычислений программа DE23 вызывает программу AIG2R обращения матрицы.

call DE26(DS, DSJ, A, B, Y, N, Eps, nWex, nErr)

DS, DSJ, A, B, N, nWex — входные параметры;
Y, Eps — входные и выходные параметры;
nErr — выходной параметр;

DS(X, Y, DY) — процедура, которая вычисляет значения производных;
DSJ(X, Y, DYJ) — процедура, которая вычисляет значения матрицы Якоби. Частная производня от правой части i-го уравнения по j-ой переменной Y(i) запоминается в элементе DYJ(i,j) ;
Real A, B — начальная и конечная точки интегрирования;
Integer N — число промежутков интегрирования от точки A до точки B . N должно быть >=1;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке A , на выходе будет содержать решение системы в точке B . Число M задает размерность решаемой системы уравнений;
Real Eps — требуемая точность интегрирования на каждом шаге; если задано значение меньшее, чем минимально-допустимое, программа увеличит его до минимально-допусимого;
Integer nWex — задаёт режим вычисления Якобиана, может принимать значения 0, 1 или 2; 0 — Якобиан вычисляется только один раз в начале программы; 1 — вычисляется 1 раз на кждом шаге перед первой итерацией; 2 — вычисляется перед каждой итерацией, т.е. от 1 до 3 раз на каждом шаге.
Integer nErr — счётчик количества точек, при которых за 3 итерации не удалось достигнуть заданной точности по какой-либо компоненте решения.

Программа DE53 решает задчу Коши для жесткой автономной системы из m обыкновенных дифференциальных уравнений первого порядка y_i‘(x) = f_i(y₁(x), y₂(x), . y_m(x)) , i=1, 2, . m методом типа Розенброка 3-го порядка аппроксимации с автоматическим выбором величины шага интегрирования. В случае неавтономной системы y’ = f(x, y(x)) введением дополнительной переменной y’_m+1 = 1 , y_m+1(x₀) = x₀ она приводится к автономному виду.

Решение системы вычисляется по формулам:

где D_n = E + α•h•f_n′ ;
E — единичная матрица; h — шаг интегрирования;
f_n′ = ∂f(y_n)/∂y — матрица Якоби системы;
α, β₂₁, β₃₁, β₃₂, p₁, p₂, p₃ — числовые коэффициенты, определяющие свойства точности и устойчивости алгоритма.

В каждом узле интегрирования программа вычисляет погрешность полученного решения. Для оценки погрешности используется идея вложенных методов, а именно, дополнительно вычисляется
y_n+1/2 = y_n + b₁k₁ + b₂k₂
Теперь оценить погрешность можно по фомуле
ε_n(2) ≈ (1/c)·D_n -1 ·(y_n+1 + y_n+1/2)
где c — коэффициент, D_n -1 — матрица, обратная к матрице D_n .

В программе вычисляется коэффициент q₂ из формулы
q₂ 3 ·║ε_n(2)║ = C·ε
где ε — заданная пользователем точность вычислений, C — коэффициент. Если q₂ , то программа повторяет вычисления из прежней точки с меньшей величиной шага. Если q₂ > 1 , то программа считает, что требуемая точность получена и переходит к вычислению следующего шага. В обоих случаях новый шаг расчитывается по формуле h_new = k·q₂·h_old , где k — понижающий коэффициент.

Норма ║ξ║ расчитывается по формуле
║ξ║ = max_1≤j≤m<|ξ_j|/(|y_n j | + r)>
где y_n j — j-ая компонента решения, r – малый положительный параметр. Если |y_n j | , применяется абсолютная погрешность, иначе относительная погрешность.

call DE53R(DSA, DSJA, X, Xout, H, Y, Eps, nWex, Error)

DSA, DSJA, Xout, nWex — входные параметры;
X, H, Y, Eps — выходные и выходные параметры;
Error — выходной параметр;

DSA(Y, DY) — процедура, которая вычисляет значения производных для автономной системы;
DSJA(Y, DYJ) — процедура, которая вычисляет значения матрицы Якоби для автономной системы. Частная производня от правой части i-го уравнения по j-ой переменной Y(i) запоминается в элементе DYJ(i,j) ;
Real X, Xout — начальная и конечная точки интегрирования, после удачного (не аварийного) завершения прграммы X = Xout ;
Real H — начальный шаг интегрирования, с которым программе предлагается начать вычисления. Допустимо задать абсолютное значение, даже если шаг должен быть отрицательным. На выходе переменная содержит величину шага в конце промежутка интегрирования;
Real Y(1:M) — массив, который на входе должен содержать начальные значения функций в точке X , на выходе будет содержать решение системы в точке Xout . Число M задает размерность решаемой системы уравнений;
Real Eps — требуемая точность интегрирования на каждом шаге. Так как алгоритм имеет 3-й порядок аппроксимации, оптимальным значением будет ε = 10 -4 ;
Integer nWex — задаёт режим вычисления Якобиана, может принимать значения 0 или 1; 0 — Якобиан вычисляется только один раз в начале программы; 1 — вычисляется на каждом шаге.
Integer Error — индикатор ошибки.
Error=0 , если интегрирование системы успешно закончено и решение вычислено с заданной точностью Eps на каждом шаге;
Error=1 , если интегрирование нельзя начать, т.к. значение |Xout-X| меньше минимально-допустимой величины шага h_min , которая вычисляется программой. Параметр H приобретает значение h_min ;
Error=2 , если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины ε_min , которая вычисляется программой. Параметр Eps приобретает значение ε_min ;
Error=65 , если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps . Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
В процессе вычислений программа DE53 вызывает программу AIG2R .

Программа DT10 решает краевую задачу для линейного дифференциального уравнения 2-го порядка y»(x) + q(x)y'(x) + p(x)y(x) = f(x) на промежутке [a, b] с граничными условиями y(a)=y_a, y(b)=y_b методом прогонки. Первая и вторая производные заменяются разностными схемами y’ = (y_i+1 — y_i-1)/2h и y» = (y_i+1 — 2y_i + y_i-1)/h² , i = 1, . n-1; x₀ = a , x_n = b ; h = x_i+1 — x_i — шаг сетки. Подстановка прводит к системе из n-2 уравнений:
a_i·y_i-1 + b_i·y_i + c_i·y_i+1 = d_i, i = 1, . n-1
где a_i = 2 — h·q(x_i) , b_i = -2(2 — h²·p(x_i)) , c_i = 2 + h·q(x_i) , d_i = 2h²·f(x_i) . Два недостающих уравнения беруться из граничных условий y(x₀) = y₀ и y(x_n) = y_n .
Получившаяся система из n уравнений с трёхдиагональной матрицей решается стандартным методом прогонки.

call DT10R(Q, P, F, A, B, YA, YB, X, Y)

Q, P, F, A, B, YA, YB — входные параметры;
X, Y — выходные параметры;

Real Q(x), P(x), F(x) — функции, определяющие дифференциальное уравнение;
Real A, B — граничные точки промежутка интегрирования;
Real YA, YB — значение фунуции в граничных точках;
Real X(1:n), Y(1:n) — массивы, которые по завершении работы программы содержат значения аргумента и функции. Число n задает разрядность сетки.

Программа DT11 решает краевую задачу для линейного дифференциального уравнения 2-го порядка y»(x) + q(x)y'(x) + p(x)y(x) = f(x) на промежутке [a, b] со смешанными граничными условиями методом прогонки. Первая и вторая производные заменяются разностными схемами y’ = (y_i+1 — y_i-1)/2h и y» = (y_i+1 — 2y_i + y_i-1)/h² , i = 1, . n-1; x₀ = a , x_n = b ; h = x_i+1 — x_i — шаг сетки. Подстановка прводит к системе из n-2 уравнений:
a_i·y_i-1 + b_i·y_i + c_i·y_i+1 = d_i, i = 1, . n-1
где a_i = 2 — h·q(x_i) , b_i = -2(2 — h²·p(x_i)) , c_i = 2 + h·q(x_i) , d_i = 2h²·f(x_i) . Два недостающих уравнения берём из смешанных граничных условий:
y'(x₀) + α·y(x₀) = α₁
y'(x_n) + β·y(x_n) = β₁
Данные граничные условия заменяем разностными схемами:
(1/2h)·(y₁ — y_-1) + α·y₀ = α₁
(1/2h)·(y_n+1 — y_n-1) + β·y_n = β₁
из которых выражаем неизвестные y_-1 и y_n+1 :
y_-1 = y₁ — 2h·(α₁ — α·y₀)
y_n+1 = y_n-1 + 2h·(β₁ — β·y_n)
Подставляем их в уравнения:
a₀·y_-1 + b₀·y₀ + c₀·y₁ = d₀
a_n·y_n-1 + b_n·y_n + c_n·y_n+1 = d_n
и получаем два недостающих уравнения для решения системы:
(b₀ + 2h·α·a₀)·y₀ + (a₀ + c₀)·y₁ = d₀ + 2h·α₁·a₀
(a_n + c_n)·y_n-1 + (b_n — 2h·β·c_n)·y_n = d_n — 2h·β₁·c_n
Получившаяся система из n уравнений с трёхдиагональной матрицей решается стандартным методом прогонки.

call DT11R(Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1, X, Y)

Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1 — входные параметры;
X, Y — выходные параметры;

Real Q(x), P(x), F(x) — функции, определяющие дифференциальное уравнение;
Real A, B — граничные точки промежутка интегрирования;
Real Alpha, Beta — коэффициент при функциях в граничных условиях;
Real Alpha1, Beta1 — правые части граничных условий;
Real X(1:n), Y(1:n) — массивы, которые по завершении работы программы содержат значения аргумента и функции. Число n задает разрядность сетки.

Видео:6. Дифференциальные уравнения, приводящиеся к однороднымСкачать

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

*****************************************************************
Описание содержимого наборов данных
ПАКЕТА НАУЧНЫХ ПРОГРАММ на языке Фортран IV
( он же пакет программ для научных исследований фирмы IBM )

Источник: ПНП Фортран IV. Общее описание.
ПРО.309.004 Д
Эстонское НПО ВТИ 1983г.
Веселовский В.Н.
*****************************************************************

PNPFS — имя папки, содержащей подпрограммы пакета в виде
исходных модулей на языке Фортран IV ДОС ЕС ЭВМ.

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

OPER.DATA — имя последовательного набора данных,
содержащего входные данные для программы-примера OPER.

BOL.DATA — имя последовательного набора данных,
содержащего входные данные для программы-примера BOL.

****************************************************************
** **
** ОГЛАВЛЕНИЕ ПАКЕТА НАУЧНЫХ ПРОГРАММ (ПНП) **
** НА ЯЗЫКЕ ФОРТРАН-4. **
** ( он же пакет программ для научных исследований фирмы IBM )**
** **
** Описание пакета имеется в книге: **
** Сборник научных программ на Фортране. Пер. с англ.- **
** М.: Статистика, 1974. **
** Примеры использования программ в расчетах приведены в **
** книге: **
** Шуп Т. Решение инженерных задач на ЭВМ: Практи- **
** ческое руководство. Пер. с англ.-М.: Мир, 1982.- **
** 238 с., ил. **
** **
****************************************************************

ЕС ЭВМ. Пакет научных программ на языке Фортран-4.
Руководство программиста. Книга 1. ПРО.309.004 Д1.

MCPY Д1.02.01.01. Размещение матрицы в памяти. Копирование мат-
рицы.
RCPY Д1.02.02.01. Размещение матрицы в памяти. Представление
строки матрицы в виде вектора.
CCPY Д1.02.03.01. Представление столбца матрицы в виде векто-
ра.
DCPY Д1.02.04.01. Представление диоганали матрицы в виде век-
тора.
XCPY Д1.02.05.01. Формирование подматрицы из данной матрицы.

MSTR Д1.02.06.01. Преобразование вида памяти по матрицы.

LOC Д1.02.07.01. Размещене матрицы в памяти. Обращене к эле-
менту массива. Определение индекса и сжатие
матрицы.
CONVT Д1.02.08.01. Преобразоване элементов матрицы из обычной
точности в удвоенную и наоборот.
ARRAY Д1.02.09.01. Преобразование одномерного массива в двумер-
ный и наоборот.

GMADD Д1.03.01.01. Сложение двух матриц общего вида.

GMSUB Д1.03.02.01. Вычитание двух матриц общего вида.

GMPRD Д1.03.03.01. Произведение двух матриц общего вида.

DGMPRD Д1.03.03.02. ( Удвоенная точность ).

GMTRA Д1.03.04.01. Транспонирование матрицы общего вида.

DGMTRA Д1.03.04.02. ( Удвоенная точность ).

GTPRD Д1.03.05.01. Умножение транспонированной матрицы на
матрицу общего вида.
DGTPRD Д1.03.05.02. ( Удвоенная точность ).

MADD Д1.03.06.01. Сложение двух матриц.

MSUB Д1.03.07.01. Вычитание двух матриц.

MPRD Д1.03.08.01. Произведение матриц (строка на столбец).

MTRA Д1.03.09.01. Транспонирование матриц.

TPRD Д1.03.10.01. Умножение транспонированной матрицы на
матрицу.
MATA Д1.03.11.01. Умножение слева матрицы на транспонирован-
ную матрицу.
SADD Д1.03.12.01. Сложение скаляра с каждым элементом матри-
цы.
SSUB Д1.03.13.01. Вычитание скаляра из каждого элемента мат
рицы.
SMPY Д1.03.14.01. Умножение матрицы на скаляр.

SDIV Д1.03.15.01. Деление матрицы на скаляр.

SCLA Д1.03.16.01. Замена каждого элемента матрицы на задан-
ный скаляр.
DCLA Д1.03.17.01. Замена диоганальных элементов матрицы на за-
данный скаляр.
RADD Д1.03.18.01. Сложение строки одной матрицы со строкой
другой матрицы.
CADD Д1.03.19.01. Сложение столбца одной матрицы со столбцом
другой матрицы.
SRMA Д1.03.20.01. Сложение одной строки матрицы с другой, ум-
ноженной на скаляр.
SCMA Д1.03.21.01. Сложение одного столбца матрицы с другим,
умноженным на скаляр.
RINT Д1.03.22.01. Перестановка двух строк матрицы.

CINT Д1.03.23.01. Перестановка двух столбцов матрицы.

RSUM Д1.03.24.01. Суммирование каждой строки матрицы.

CSUM Д1.03.25.01. Суммирование каждого столбца матрицы.

RTAB Д1.03.26.01. Табулирование строк матрицы.

CTAB Д1.03.27.01. Табулирование столбцов матрицы.

RSRT Д1.03.28.01. Сортировка строк матрицы.

CSRT Д1.03.29.01. Сортировка столбцов матрицы.

RCUT Д1.03.30.01. Разбиение матрицы на две по строке.

CCUT Д1.03.31.01. Разбиение матрицы на две по столбцу.

RTIE Д1.03.32.01. Объединение двух матриц по строке.

CTIE Д1.03.33.01. Объединение двух матриц по столбцу.

MPRC Д1.03.34.01. Перестановка строк или столбцов матрицы.

DMPRC Д1.03.34.02. ( Удвоенная точность ).

MFUN Д1.03.35.01. Преобразоване матрицы с помощью функции.

RECP Д1.03.36.01. Вычисление обратной величины элемента.

Д1.04 Обращене матриц, системы алгебраи-
ческих уравнений и родственные темы.
———————————————

MINV Д1.04.01.01. Обращение матрицы.

DMINV Д1.04.01.02. ( Удвоенная точность ).

SINV Д1.04.02.01. Обращение симметричной положительно опре-
деленной матрицы.
DSINV Д1.04.02.02. ( Удвоенная точность ).

SIMQ Д1.04.03.01. Решение системы линейных алгебраических
уравнений методом мсключения.
DSIMQ Д1.04.03.02. ( Удвоенная точность ).

GELG Д1.04.04.01. Решение системы линейных уравнений общего
вида с правыми частями методом исключения
DGELG Д1.04.04.02. ( Удвоенная точность ). Гаусса.

RSLMC Д1.04.05.01. Решение системы линейных уравнений с ите-
рационным уточнением.
FACTR Д1.04.06.01. Разложение вырожденной матрицы на произведе-
ние двух треугольных.
MFGR Д1.04.07.01. Разложение матрицы на множители и опреде-
ление ранга.
DMFGR Д1.04.07.02 ( Удвоенная точность ).

GELS Д1.04.08.01. Решение системы линейных уравнений с сим-
метричной матрицей коэффициентов.
DGELS Д1.04.08.02. ( Удвоенная точность ).

GELB Д1.04.09.01. Решение системы линейных уравнений с ленточ-
ной матрицей коэффициентов.
DGELB Д1.04.09.02. ( Удвоенная точность ).

MTDS Д1.04.10.01. Деление матрицы на треугольную матрицу.

DMTDS Д1.04.10.02. ( Удвоенная точность ).

MLSS Д1.04.11.01. Решение системы линейных уравнений с сим-
метричной положительной полуопределенной
DMLSS Д1.04.11.02. ( Удвоенная точность ). матрицей.

MCHB Д1.04.12.01. Разложение на множители симметричной поло-
жительно определенной матрицы.
DMCHB Д1.04.12.02. ( Удвоенная точность ).

MFSS Д1.04.13.01. Разложение на множители и определение ранга
симметричной положительно опрнделенной
DMFSS Д1.04.13.02. ( Удвоенная точность ). матрицы.

MFSD Д1.04.14.01. Разложение на множители симметричной поло-
житнльно определенной матрицы.
DMFSD Д1.04.14.02. ( Удвоенная точность ).

LLSQ Д1.04.15.01. Решени линейной задачи методом наименьших
квадратов.
DLLSQ Д1.04.15.02. ( Удвоенная точность ).

Д1.05 Собственные значения, собственные векторы и
родственные темы.
—————————————————

EIGEN Д1.05.01.01. Вычисление собственных значений и собствен-
ных векторов действительной симметричной
матрицы.
DEIGEN Д1.05.01.02. ( Удвоенная точность ).

NROOT Д1.05.02.01. Вычисление собственных значений и собствен-
ных векторов действительной несимметричной
матрицы вида (В)**-1 * А.
DNROOT Д1.05.02.02. ( Удвоенная точность ).

ATEIG Д1.05.03.01. Вычисление собственных значений действитель-
ной почти треугольной матрицы.
HSBG Д1.05.04.01. Приведение действительной матрицы к верхней
почти треугольной форме.

MATIN Д1.06.01.01. Чтение с перфокарт в память матрицы общего
вида, симметричной или диоганальной.
DMATIN Д1.06.01.02. ( Удвоенная точность ).

MXOUT Д1.06.02.01. Печать матрицы общего вида, симметричной
или диоганальной.
DMXOUT Д1.06.02.02. ( Удвоенная точность ).

ЕС ЭВМ. Пакет научных программ на языке Фортран-4.
Руководство программиста. Книга 2. ПРО.309.004 Д2.

ABSNT Д2.01.01.01. Выявление недостающих данных ( 0 элементы ).

SUBST Д2.01.02.01. Выбор подмножества из матрицы наблюдений.

TALLY Д2.01.03.01. Числовые характеристики статистического ряда.

BOUND Д2.01.04.01. Отбор наблюдений в пределах заданных границ.

TAB1 Д2.01.05.01. Вычисление частот попадания переменной в за-
данный интервал.
TAB2 Д2.01.06.01. Вычисление частот попадания для двух перемен-
ных.
SUBMX Д2.01.07.01. Построение матрицы подмножества по условиям.

CORRE Д2.02.01.01. Вычисление коэффициента корреляции.

MISR Д2.02.02.01. Вычисление коэффициентов корреляции и линии
регресси.
ORDER Д2.02.03.01. Образоване подмножества из большой корреля-
ционной матрицы.
MULTR Д2.02.04.01. Множественный линейный регрессионный анализ
переменных.
GDATA Д2.02.05.01. Образование матрицы данных для полиноминой
регресси.
STPRG Д2.02.06.01. Множественная линейная регрессия с выбором
существенных факторов.
PROBT Д2.02.07.01. Пробит — анализ.

CANOR Д2.02.08.01. Каноническая корреляция между двумя мно-
жествами переменных.

AVDAT Д2.03.01.01. Размещение данных для дисперсного анализа.

AVCAL Д2.03.02.01. Вычисление операторов факторного экспери-
мента.
MEANQ Д2.03.03.01. Вычисление сумм квадратов отклонений.

Д2.04 Дискриминантный анализ.

DMATX Д2.04.01.01. Вычисление средних и ковариационной матрицы.

DISCR Д2.04.02.01. Вычисление дискриминантных функций.

TRACE Д2.05.01.01. Вычисление накопленных отклонений собствен-
ных значений корреляционной матрицы.
LOAD Д2.05.02.01. Выччисление матрицы факторных нагрузок

VARMX Д2.05.03.01. Ортогональное вращение матрицы факторов.

AUTO Д2.06.01.01. Нахождение автоковариаций ряда А для за-
паздывания аргумента от 0 до L-1.
CROSS Д2.06.02.01. Вычисление взаимоковариационной функции.

SMO Д2.06.03.01. Взвешивание мсходных данных.

EXSMO Д2.06.04.01. Тройное экспоненциальное сглаживание.

KOLMO Д2.07.01.01. Проверка разностимежду теоретическим и эм-
пирическим распределениями с использованием
критерия Колмогорова-Смирнова.
KOLM2 Д2.07.02.01. Критерий Кломогорова-Смирнова для двух вы-
борок.
SMIRN Д2.07.03.01. Вычисление значения предела функции распре-
деления Колмогорова-Смирнова.
CHISQ Д2.07.04.01. Вычисление критерия для таблицы сопряжен-
ности признаков.
KRANK Д2.07.05.01. Вычисление коэффициента ранговой корреляции
по Кендаллу.
MPAIR Д2.07.06.01. Критерий Вилкоксона.

QTEST Д2.07.07.01. Отличие групп по критерию Q Кохрена.

RANK Д2.07.08.01. Расположение элементов вектора в возраста-
ющем порядке.
SIGNT Д2.07.09.01. Выполняет непараметрический критерий знаков.

SRANK Д2.07.10.01. Вычисление коэффициента ранговой корреляции
по Смирмену.
TIE Д2.07.11.01. Вычисление поправочного коэффициента для
ранговой корреляции.
TWOAV Д2.07.12.01. Двухфакторный дисперсный анализ по Фридману.

UTEST Д2.07.13.01. Критерий Минна-Уитни.

WTEST Д2.07.14.01. Коэффициент соотвенствия Кендалла.

RANDU Д2.08.01.01. Вычисление равномерно распределенног случай-
ного числа.
GAUSS Д2.08.02.01. Вычисление нормально распределенного случай-
ного числа.
NDTR Д2.08.03.01. Вычисление вероятности.

BDTR Д2.08.04.01. Бета-распределение.

CDTR Д2.08.05.01. Хи квадрат-распределение.

NDTRI Д2.08.06.01. Вычисление значения аргумента по заданному
значению вероятности.

MOMEN Д2.09.01.01. Вычисление первых четырех моментов.

TTEST Д2.09.02.01. Вычисление t-статистики.

BISER Д2.09.03.01. Вычисление бисериального коэффициента
корреляции.
PHI Д2.09.04.01. Вычисление Ф-коэффициента.

POINT Д2.09.05.01. Вычислене точечно-бисериального коэффи-
циента корреляции.
TETRA Д2.09.06.01. Вычисление тетрахорического коэффициента
корреляции.
SRATE Д2.09.07.01. Коэффициент выживания.

RK1 Д2.10.01.01. Решение обыкновенного дифференциального урав-
нения первого порядка методом Рунге-Кутта с
выдачей решения в конечной точке.
RK2 Д2.10.02.01. Решение обыкновенного дифференциального урав-
нения первого порядка методом Рунге-Кутта с
выдачей таблицы значений решения.
RKGS Д2.10.03.01. Решение системы обыкновенных дифференциаль-
ных уравнений первого порядка методом Рун-
ге-Кутта.
DRKGS Д2.10.03.02. ( Удвоенная точность ).

HPCG Д2.10.04.01. Решение системы обыкновенных дифференциаль-
ных уравнений первого порядка методом Хем-
минга.
DHPCG Д2.10.04.02. ( Удвоенная точность ).

HPCL Д2.10.05.01. Решение системы линейных дифференциальных
уравнений первого порядка методом Хемминга.
DHPCL Д2.10.05.02. ( Удвоенная точность ).

LBVP Д2.10.06.01. Решение системы линейных обыкновенных диф-
ференциальных уравнений первого порядка с
линейными краевыми условиями.
DLBVP Д2.10.06.02. ( Удвоенная точность ).

ЕС ЭВМ. Пакет научных программ на языке Фортран-4.
Руководство программиста. Книга 3. ПРО.309.004 Д3.

PADD Д3.01.01.01. Сложение двух полиномов.

DPADD Д3.01.01.02. ( Удвоенная точность ).

PSUB Д3.01.02.01. Вычитание одного полинома из другого.

DPSUB Д3.01.02.02. ( Удвоенная точность ).

PMPY Д3.01.03.01. Умножение двух полиномов.

DPMPY Д3.01.03.02. ( Удвоенная точность ).

PDIV Д3.01.04.01. Деление одного полинома на другой.

DPDIV Д3.01.04.02. ( Удвоенная точность ).

PCLA Д3.01.05.01. Замена одного полинома другим.

DPCLA Д3.01.05.02. ( Удвоенная точность ).

PADDM Д3.01.06.01. Сложение одного полинома с другим, умно-
женным на число.
DPADDM Д3.01.06.02. ( Удвоенная точность ).

PVAL Д3.01.07.01. Вычисление значения полинома.

DPVAL Д3.01.07.02. ( Удвоенная точность ).

PVSUB Д3.01.08.01. Замена переменной полинома другим полиномом.

DPVSUB Д3.01.08.02. ( Удвоенная точность ).

PILD Д3.01.09.01. Вычисление значения полинома и его первой
производной.
DPILD Д3.01.09.02. ( Удвоенная точность ).

PDER Д3.01.10.01. Вычисление производной полинома.

DPDER Д3.01.10.02. ( Удвоенная точность ).

PINT Д3.01.11.01. Вычисление неопределенного интеграла от по-
линома.
DPINT Д3.01.11.02. ( Удвоенная точность ).

PQSD Д3.01.12.01. Деление полинома на квадратный трехчлен.

DPQSD Д3.01.12.02. ( Удвоенная точность ).

PCLD Д3.01.13.01. Вычисление коэффициентов Р(x-u) по заданным
коэффициентам P(x) и U.
DPCLD Д3.01.13.02. ( Удвоенная точность ).

PGCD Д3.01.14.01. Вычисление наибольшего общего делителя двух
полиномов.
DPGCD Д3.01.14.02. ( Удвоенная точность ).

PNORM Д3.01.15.01. Нормировка вектора коэффициентов полинома.

DPNORM Д3.01.15.02. ( Удвоенная точность ).

PECN Д3.01.16.01. Экономизация полинома для несимметричного
интервала.
DPECN Д3.01.16.02. ( Удвоенная точность ).

PECS Д3.01.17.01. Экономизация полинома для несимметричного
интервала.
DPECS Д3.01.17.02. ( Удвоенная точность ).

POLRT Д3.02.01.01. Вычисление всех корней полинома с действи-
тельными коэффициентами методом Ньютона-
Рафсона.
DPOLRT Д3.02.01.02. ( Удвоенная точность ).

PRQD Д3.02.02.01. Вычисление корней полинома с действительными
коэффициентами посредством QD-алгоритма.
DPRQD Д3.02.02.02. ( Удвоенная точность ).

PRBM Д3.02.03.01. Вычисление корней полинома с действительными
коэффициентами методом Берстоу.
DPRBM Д3.02.03.02. ( Удвоенная точность ).

PQFB Д3.02.04.01. Определение квадратного множителя полинома
с действительными коэффициентами.
DPQFB Д3.02.04.02. ( Удвоенная точность ).

CNP Д3.03.01.01. Вычисление значений полиномов Чебышева.

DCNP Д3.03.01.02. ( Удвоенная точность ).

CNPS Д3.03.02.01. Вычисление значений функции, разложенной в
ряд по полиномам Чебышева.
DCNPS Д3.03.02.02. ( Удвоенная точность ).

TCNP Д3.03.03.01. Преобразование разложения функции по полино-
мам Чебышева в полином.
DTCNP Д3.03.03.02. ( Удвоенная точность ).

CSP Д3.03.04.01. Вычисление значений смещенных полиномов
Чебышева.
DCSP Д3.03.04.02. ( Удвоенная точность ).

CSPS Д3.03.05.01. Вычисление функции, разложенной в ряд по
смещенным полиномам Чебышева.
DCSPS Д3.03.05.02. ( Удвоенная точность ).

TCSP Д3.03.06.01. Преобразование разложения функции по смещен-
ным полиномам Чебышева в полином.
DTCSP Д3.03.06.02. ( Удвоенная точность ).

HEP Д3.03.07.01. Вычисление значений полиномов Эрмита.

DHEP Д3.03.07.02. ( Удвоенная точность ).

HEPS Д3.03.08.01. Вычисление значения функции, разложенной в
ряд по полиномам Эрмита.
DHEPS Д3.03.08.02. ( Удвоенная точность ).

THEP Д3.03.09.01. Преобразование разложения функции по поли-
номам Эрмита в полином.
DTHEP Д3.03.09.02. ( Удвоенная точность ).

LAP Д3.03.10.01. Вычисление значений полиномов Лагерра.

DLAP Д3.03.10.02. ( Удвоенная точность ).

LAPS Д3.03.11.01. Вычисление значения функции, разложенной в
ряд по полиномам Лагерра.
DLAPS Д3.03.11.02. ( Удвоенная точность ).

TLAP Д3.03.12.01. Преобразование разложения функции по поли-
номам Лагерра в полином.
DTLAP Д3.03.12.02. ( Удвоенная точность ).

LEP Д3.03.13.01. Вычислениние значений полиномов Лежандра.

DLEP Д3.03.13.02. ( Удвоенная точность ).

LEPS Д3.03.14.01. Вычисление значений функции, разложенной в
ряд по полиномам Лежандра.
DLEPS Д3.03.14.02 ( Удвоенная точность ).

TLEP Д3.03.15.01. Преобразование разложения функции по поли-
номом Лежандра в полином.
DTLEP Д3.03.15.02. ( Удвоенная точность ).

RTWI Д3.04.01.01. Уточнение корня нелинейного уравнения ите-
рационным методом Вегстейна.
DRTWI Д3.04.01.02. ( Удвоенная точность ).

RTMI Д3.04.02.01. Вычисление корня нелинейного уравнения внут-
ри интервала итерационным методом Мюллера.
DRTMI Д3.04.02.02. ( Удвоенная точность ).

RTNI Д3.04.03.01. Уточнение корня нелинейного уравнения мето-
дом Ньютона.
DRTNI Д3.04.03.02. ( Удвоенная точность ).

FMFP Д3.05.01.01. Отыскание локального минимума функции нес-
кольких переменных методом Флетчера-Павелла.
DFMFP Д3.05.01.02. ( Удвоенная точность ).

FMCG Д3.05.02.01. Отыскание локального минимума функции нес-
кольких переменных методом сопряженных гра-
диентов.
DFMCG Д3.05.02.02. ( Удвоенная точность ).

PPRCN Д3.06.01.01. По данным векторам подстановок вычисление
произведения и сопряжения.
PERM Д3.06.02.01. Операции с подстановками и транспозициями.

TEAS Д3.07.01.01. Вычисление предела данной последователь-
ности с помощью эпсилон алгоритма.
DTEAS Д3.07.01.02. ( Удвоенная точность ).

TEUL Д3.07.02.01. Вычисление суммы функциональной последова-
тельности.
DTEUL Д3.07.02.02. ( Удвоенная точность ).

GMMMA Д3.08.01.01. Вычисление гамма-функции.

DLGAM Д3.08.02.01. Вычисление натурального логарифма гамма-
функции ( Удвоенная точность ).
BESJ Д3.08.03.01. Вычисление функции Бесселя Jn(X).

BESY Д3.08.04.01. Вычисление функции Бесселя Yn(X).

BESI Д3.08.05.01. Вычисление функции Бесселя In(X).

BESK Д3.08.06.01. Вычисление функции Бесселя Kn(X).

EXPI Д3.08.07.01. Вычисление показательного интеграла.

SICI Д3.08.08.01. Вычисление интегрального синуса и косинуса.

CS Д3.08.09.01. Вычисление интегралов Френеля.

CEL1 Д3.08.10.01. Вычисление полного эллиптического интеграла
первого рода.
DCEL1 Д3.08.10.02. ( Удвоенная точность ).

CEL2 Д3.08.11.01. Вычисление полного эллиптического интеграла
второго рода.
DCEL2 Д3.08.11.02. ( Удвоенная точность ).

ELI1 Д3.08.12.01. Вычисление эллиптического интеграла первого
рода.
DELI1 Д3.08.12.02. ( Удвоенная точность ).

ELI2 Д3.08.13.01. Вычисление обобщеного эллиптического интег-
рала второго рода.
DELI2 Д3.08.13.02. ( Удвоенная точность ).

JELF Д3.08.14.01. Вычисление трех эллиптических функций Якоби
SN, CN и DN.
DJELF Д3.08.14.02. ( Удвоенная точность ).

ЕС ЭВМ. Пакет научных программ на языке Фортран-4.
Руководство программиста. Книга 4. ПРО.309.004 Д4.

DGT3 Д4.01.01.01. Дифференцирование функции, заданной таблицей
значений в неравноотстоящих точках
DDGT3 Д4.01.01.02. ( Удвоенная точность ).

DET3 Д4.01.02.01. Дифференцирование функции, заданной таблицей
значений в равноотстоящих точках, по формуле
Лежандра с тремя узлами.
DDET3 Д4.01.02.02. ( Удвоенная точность ).

DET5 Д4.01.03.01. Дифференцирование функции, заданной таблицей
значений в равноотстоящих точках, по формуле
Лежандра с пятью узлами.
DDET5 Д4.01.03.02. ( Удвоенная точность ).

DCAR Д4.01.04.01. Дифференцирование функции в центре интервала
методом экстрополяции Ричардсона и Ромберга.
DDCAR Д4.01.04.02. ( Удвоенная точность ).

DBAR Д4.01.05.01. Дифференцирование функции на границе интер-
вала методом экстрополяции Ричардсона и
Ромберга.
DDBAR Д4.01.05.02. ( Удвоенная точность ).

QTFG Д4.02.01.01. Интегрирование функции, заданной таблицей
значений в неравноотстоящих точках, по прави-
лу трапеций.
DQTFG Д4.02.01.02. ( Удвоенная точность ).

QTFE Д4.02.02.01. Интегрирование функции, заданной таблицей
значений в равноотстоящих точках, по прави-
лу трапеций.
DQTFE Д4.02.02.02. ( Удвоенная точность ).

QSF Д4.02.03.01. Интегрирование функции, заданной таблицей
значений в равноотстоящих точках, по прави-
лу Ньютона-Котеса.
DQSF Д4.02.03.02. ( Удвоенная точность ).

QHFG Д4.02.04.01. Интегрирование функции, заданной таблицей
значений функции и ее первой производной в
неравноотстоящих точках, по правилу Эрмита
первого порядка.
DQHFG Д4.02.04.02. ( Удвоенная точность ).

QHFE Д4.02.05.01. Интегрирование функции, заданной таблицей
значений функции и ее первой производной в
равноотстоящих точках, по правилу Эрмита
первого порядка.
DQHFE Д4.02.05.02. ( Удвоенная точность ).

QHSG Д4.02.06.01. Интегрирование функции, заданной таблицей
значений функции и ее первой и второй произ-
водными в неравноотстоящих точках, по правилу
Эрмита второго порядка.
DQHSG Д4.02.06.02. ( Удвоенная точность ).

QHSE Д4.02.07.01. Интегрирование функции, заданной таблицей
значений функции и ее первой и второй произ-
водными в равноотстоящих точках, по правилу
Эрмита второго порядка.
DQHSE Д4.02.07.02. ( Удвоенная точность ).

QATR Д4.02.08.01. Интегрирование функции по правилу трапеций
с экстрополяцией по методу Ромберга.
DQATR Д4.02.08.02. ( Удвоенная точность ).

QG2 Д4.02.09.01.=:
QG3 Д4.02.09.02. :
QG4 Д4.02.09.03. :
QG5 Д4.02.09.04. : Интегрирование функции
QG6 Д4.02.09.05. :==== по квадратурной формуле Гаусса
QG7 Д4.02.09.06. : Гаусса с N узлами.
QG8 Д4.02.09.07. :
QG9 Д4.02.09.08. :
QG10 Д4.02.09.09.=:

DQG4 Д4.02.09.10.=:
DQG8 Д4.02.09.11. :
DQG12 Д4.02.09.12. :=== ( Удвоенная точность ).
DQG16 Д4.02.09.13. :
DQG24 Д4.02.09.14. :
DQG32 Д4.02.09.15.=:

QL2 Д4.02.10.01.=:
QL3 Д4.02.10.02. :
QL4 Д4.02.10.03. : Вычисление несобсвенного интеграла функции
QL5 Д4.02.10.04. : 0 -x
QL6 Д4.02.10.05. :=== _/- e f(x) dx по квадратурной формуле
QL7 Д4.02.10.06. : бескон.
QL8 Д4.02.10.07. : Гаусса-Лагерра с N узлами.
QL9 Д4.02.10.08. :
QL10 Д4.02.10.09.=:

DQL4 Д4.02.10.10.=:
DQL8 Д4.02.10.11. :
DQL12 Д4.02.10.12. :=== (Удвоенная точность ).
DQL16 Д4.02.10.13. :
DQL24 Д4.02.10.14. :
DQL32 Д4.02.10.15.=:

QH2 Д4.02.11.01.=:
QH3 Д4.02.11.02. : Вычисление несобственного интеграла
QH4 Д4.02.11.03. : +бескон. -x2
QH5 Д4.02.11.04. :=== y= / e f(x) dx по
QH6 Д4.02.11.05. : -бескон.
QH7 Д4.02.11.06. :
QH8 Д4.02.11.07. : квадратурной формуле Гаусса-Эрмита с
QH9 Д4.02.11.08. : N узлами.
QH10 Д4.02.11.09.=:

DQH8 Д4.02.11.10.=:
DQH16 Д4.02.11.11. :
DQH24 Д4.02.11.12. :=== ( Удвоенная точность ).
DQH32 Д4.02.11.13. :
DQH48 Д4.02.11.14. :
DQH64 Д4.02.11.15.=:

QA2 Д4.02.12.01.=:
QA3 Д4.02.12.02. : Вычисление несобственного интеграла
QA4 Д4.02.12.03. : бескон. -x
QA5 Д4.02.12.04. :=== y= / (e f(x) / sqrt(x) ) dx по
QA6 Д4.02.12.05. : 0
QA7 Д4.02.12.06. : присоединенной квадратурной формуле
QA8 Д4.02.12.07. : Гаусса-Лагерра с N узлами.
QA9 Д4.02.12.08.=:

QA10 Д4.02.12.09.=:
DQA4 Д4.02.12.10. :
DQA8 Д4.02.12.11 :
DQA12 Д4.02.12.12. :=== ( Удвоенная точность ).
DQA16 Д4.02.12.13. :
DQA24 Д4.02.12.14. :
DQA32 Д4.02.12.15.=:

ALI Д4.03.01.01. Интерполирование значения функции для задан-
ного значения аргумента по таблице методом
Эйткина-Лагранжа.
DALI Д4.03.01.02. ( Удвоенная точность ).

AHI Д4.03.02.01. Интерполирование функции с помощью процесса
Эйткина-Эрмита.
DAHI Д4.03.02.02. ( Удвоенная точность ).

ACFI Д4.03.03.01. Интерполирование функции с помощью непрерыв-
ной дроби.
DACFI Д4.03.03.02. ( Удвоенная точность ).

ATSG Д4.03.04.01. Выборка таблицы из таблицы общего вида.

DATSG Д4.03.04.02. ( Удвоенная точность ).

ATSM Д4.03.05.01. Выборка таблицы из таблицы с монотонными аргу-
ментами.
DATSM Д4.03.05.02. ( Удвоенная точность ).

ATSE Д4.03.06.01. Выборка таблицы из таблицы с равноотстоящими
аргументами.
DATSE Д4.03.06.02. ( Удвоенная точность ).

SG13 Д4.04.01.01. Сглаживание функции заданной таблицей с по-
мощью многочлена 1-ой степени по методу на-
именьших квадратов.
DSG13 Д4.04.01.02. ( Удвоенная точность ).

SE13 Д4.04.02.01. Сглаживание функции, заданной таблицей в
равноотстоящих точках, многочленом 1-ой сте-
пени по методу наименьших квадратов.
DSE13 Д4.04.02.02. ( Удвоенная точность ).

SE15 Д4.04.03.01. Сглаживание функции, заданной таблицей в
равноотстоящих точках, многочленом 1-ой сте-
пени по 5-ти точкам методом наименьших квад-
ратов.
DSE15 Д4.04.03.02. ( Удвоенная точность ).

SE35 Д4.04.04.01. Сглаживание функции, заданной таблицей в
равноотстоящих точках многочленом 3-ей
степени по 5-ти точкам методом наименьших
квадратов.
DSE35 Д4.04.04.02. ( Удвоенная точность ).

APFS Д4.05.01.01. Решение системы нормальных уравнений, воз-
никающих при приближении по методу наимень-
ших квадратов.
DAPFS Д4.05.01.02. ( Удвоенная точность ).

APCH Д4.05.02.01. Составление системы нормальных уравнений при
приближении по методу наименьших квадратов.
DAPCH Д4.05.02.02. ( Удвоенная точность ).

ARAT Д4.05.03.01. Аппроксимация дискретной функции рациональ-
ной функцией методом наименьших квадратов.
DARAT Д4.05.03.03. ( Удвоенная точность ).

FRAT Д4.05.03.02. Обработка данных в методе наименьших ква-
дратов.
DFRAT Д4.05.03.04. ( Удвоенная точность ).

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

APMM Д4.05.05.01. Аппроксимация дискретной функции линейной
комбинацией функций по методу Чебышева.
DAPMM Д4.05.05.02. ( Удвоенная точность ).

FORIF Д4.06.01.01. Вычисление коэффициентов Фурье периодичес-
кой функции.
FORIT Д4.06.02.01. Вычисление коэффициентов Фурье периодичес-
кой функции, заданной таблично.
HARM Д4.06.03.01. Комплексный трехмерный анализ Фурье.

DHARM Д4.06.03.02. ( Удвоенная точность ).

RHARM Д4.06.04.01. Одномерный действительный анализ Фурье.

DRHARM Д4.06.04.02. ( Удвоенная точность ).

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

Численное решение математических моделей объектов заданных системами дифференциальных уравнений

Введение:

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

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

Возникает необходимость использовать численные методы, наиболее известным из которых является метод Рунге — Кутты [1]. Что касается Python, то в публикациях по численным методам, например [2,3], данных по применение Рунге — Кутты крайне мало, а по его модификации — методу Рунге-Кутта-Фельберга вообще нет.

В настоящее время, благодаря простому интерфейсу, наибольшее распространение в Python имеет функцию odeint из модуля scipy.integrate. Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие.

Целью настоящей публикации является сравнительный анализ перечисленных средств численного решения систем дифференциальных уравнений с модифицированным автором под Python методом Рунге-Кутта-Фельберга. В публикации так же приведены решения по краевым задачам для систем дифференциальных уравнений (СДУ).

Краткие теоретические и фактические данные по рассматриваемым методам и программным средствам для численного решения СДУ

Для одного дифференциального уравнения n – го порядка, задача Коши состоит в нахождении функции, удовлетворяющей равенству:

и начальным условиям

Перед решением эта задача должна быть переписана в виде следующей СДУ

(1)

с начальными условиями

Модуль имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (задача Коши). Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач.

Функция odeint() имеет три обязательных аргумента и много опций. Она имеет следующий формат odeint(func, y0, t[,args=(), . ]) Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2. yn], а второй – имя независимой переменной.

Функция func должна возвращать список из n значений функций при заданном значении независимого аргумента t. Фактически функция func(y,t) реализует вычисление правых частей системы (1).

Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений при t=t0.

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

Функция odeint() возвращает массив размера len(t) x len(y0). Функция odeint() имеет много опций, управляющих ее работой. Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений ei для каждого значения yi по формуле

Они могут быть векторами или скалярами. По умолчанию

Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Она создает объект ОДУ (тип scipy.integrate._ode.ode). Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы. Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей.

Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций. Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши.

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

При численном решении задачи Коши

(2)

(3)

по известному решению в точке t =0 необходимо найти из уравнения (3) решение при других t. При численном решении задачи (2),(3) будем использовать равномерную, для простоты, сетку по переменной t с шагом т > 0.

Приближенное решение задачи (2), (3) в точке обозначим . Метод сходится в точке если при . Метод имеет р-й порядок точности, если , р > 0 при . Простейшая разностная схема для приближенного решения задачи (2),(3) есть

(4)

При имеем явный метод и в этом случае разностная схема аппроксимирует уравнение (2) с первым порядком. Симметричная схема в (4) имеет второй порядок аппроксимации. Эта схема относится к классу неявных — для определения приближенного решения на новом слое необходимо решать нелинейную задачу.

Явные схемы второго и более высокого порядка аппроксимации удобно строить, ориентируясь на метод предиктор-корректор. На этапе предиктора (предсказания) используется явная схема

(5)

а на этапе корректора (уточнения) — схема

В одношаговых методах Рунге—Кутта идеи предиктора-корректора реализуются наиболее полно. Этот метод записывается в общем виде:

(6),

Формула (6) основана на s вычислениях функции f и называется s-стадийной. Если при имеем явный метод Рунге—Кутта. Если при j>1 и то определяется неявно из уравнения:

(7)

О таком методе Рунге—Кутта говорят как о диагонально-неявном. Параметры определяют вариант метода Рунге—Кутта. Используется следующее представление метода (таблица Бутчера)

Одним из наиболее распространенных является явный метод Рунге—Кутта четвертого порядка

(8)

Метод Рунге—Кутта— Фельберга

Привожу значение расчётных коэффициентов метода

(9)

С учётом(9) общее решение имеет вид:

(10)

Это решение обеспечивает пятый порядок точности, остаётся его адаптировать к Python.

Вычислительный эксперимент по определению абсолютной погрешности численного решения нелинейного дифференциального уравнения с использованием обеих функций def odein(),def oden() модуля scipy.integrate и адаптированного к Python методов Рунге—Кутта и Рунге—Кутта— Фельберга

Адаптированные к Python методы Рунге—Кутта и Рунге—Кутта— Фельберга имеют меньшую абсолютную, чем решение с применением функции odeint, но большую, чем с использованием функции edu. Необходимо провести исследование быстродействия.

Численный эксперимент по сравнению быстродействия численного решения СДУ при использовании функции ode с атрибутом dopri5 (метод Рунге – Кутты 5 порядка) и с использованием адаптированного к Python метода Рунге—Кутта— Фельберга

Сравнительный анализ проведём на примере модельной задачи, приведенной в [2]. Чтобы не повторять источник, приведу постановку и решение модельной задачи из [2].

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

где – радиус вектор движущегося тела, – вектор скорости тела, – коэффициент сопротивления, вектор силы веса тела массы m, g – ускорение свободного падения.

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

К системе следует добавить начальные условия: (h начальная высота), . Положим . Тогда соответствующая система ОДУ 1 – го порядка примет вид:

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

Flight time = 1.2316 Distance = 5.9829 Height =1.8542
Flight time = 1.1016 Distance = 4.3830 Height =1.5088
Flight time = 1.0197 Distance = 3.5265 Height =1.2912
Flight time = 0.9068 Distance = 2.5842 Height =1.0240
Время на модельную задачу: 0.454787

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

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6

Функция increment(f, t, y, tau) обеспечивает пятый порядок численного метода решения. Остальные особенности программы можно посмотреть в следующем листинге:

Время на модельную задачу: 0.259927

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

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

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

(11)

Для решения задачи (11) используем следующий алгоритм:

1. Решаем первые три неоднородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

2. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

3. Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:

4. Общее решение краевой задачи (11) при помощи решений задач Коши записывается в виде линейной комбинации решений:

где p2, p3 — некоторые неизвестные параметры.

5. Для определения параметров p2, p3, используем краевые условия последних двух уравнений (11), то есть условия при x = b. Подставляя, получим систему линейных уравнений относительно неизвестных p2, p3:
(12)
Решая (12), получим соотношения для p2, p3.

По приведенному алгоритму с применением метода Рунге—Кутта—Фельберга получим следующую программу:

y0[0]= 0.0
y1[0]= 1.0
y2[0]= 0.7156448588231397
y3[0]= 1.324566562303714
y0[N-1]= 0.9900000000000007
y1[N-1]= 0.1747719838716767
y2[N-1]= 0.8
y3[N-1]= 0.5000000000000001
Время на модельную задачу: 0.070878

Вывод

Разработанная мною программа отличается от приведенной в [3] меньшей погрешностью, что подтверждает приведенный в начале статьи сравнительный анализ функции odeint с реализованным на Python метода Рунге—Кутта—Фельберга.

3. Н.М. Полякова, Е.В. Ширяева Python 3. Создание графического интерфейса пользователя (на примере решения методом пристрелки краевой задачи для линейных обыкновенных дифференциальных уравнений). Ростов-на-Дону 2017.

🎬 Видео

Решение дифференциальных уравнений. Практическая часть. 11 класс.Скачать

Решение дифференциальных уравнений. Практическая часть. 11 класс.Скачать

Решение дифференциальных уравнений. Практическая часть. 11 класс.Скачать

Решить дифференциальное уравнение #laplacetransform #differentialequation #equation #derivativesСкачать

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

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

13. Как решить дифференциальное уравнение первого порядка?Скачать

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

Линейное неоднородное дифференциальное уравнение с постоянными коэффициентами 4y''-y=x^3-24x #1Скачать

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

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

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

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