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

Обыкновенные дифференциальные уравнения

Программа 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 задает разрядность сетки.

О.Арушанян Численные решения обыкновенных дифференциальных уравнений на Фортране

1

Подробное описание

Арушанян О.Б., Залеткин С.Ф. Численное решение обыкновенных дифференциальных уравнений на Фортране. М. МГУ 1990г. 336 с. Твердый переплет, 15*22 см. (ISBN: 5-211-00957-6 / 5211009576) В книге рассматриваются способы решения задачи Коши для уравнений первого и второго порядков конечно-разностными методами, позволяющие оценивать погрешность вычислений и произвольно изменять шаг интегрирования, эффективные матричные методы решения линейных систем с постоянными коэффициентами, методы решения жестких задач. Приводятся тексты программ на языке Фортран, реализующих различные численные алгоритмы и составляющих функциональное наполнение пакета программ для решения обыкновенных дифференциальных уравнений, разработанного в НИВЦ МГУ. Для инженеров и научных работников, чья работа связана с решением дифференциальных уравнений. Тираж книги 9225 экз. Состояние: отличное Минимальная сумма покупки 150 рублей . Пишу для самых «продвинутых» покупателей : ЕСЛИ книга стоит меньше 150 р. и Вы делаете заказ, то она будет стоить 150 р. . ЗАДАВАЙТЕ ВСЕ ВОПРОСЫ ДО ТОГО КАК СДЕЛАЛИ ЗАКАЗ . Уважаемые покупатели после того как Вы сделали денежный перевод, сообщайте об этом на сайте . Если для Вас важно наличие-отсутствие библиотечного штампа, задавайте об этом вопрос . ВНИМАНИЕ! Просьба уточнять о наличии лота до заказа. В страны бывшие СНГ не отправляю, кроме прибалтики . Щеглов Максим Игоревич Факультет: Вычислительной техники и информатики Специальность: Программное обеспечение автоматизированных систем Тема магистерской диссертации: «Анализ и оценка эффективности параллельных многошаговых блочных методов решения ОДУ на кластере» Научный руководитель: Фельдман Лев Петрович профессор, д.т.н. Библиотека Авторы: Л.П. Фельдман, И.А. Назарова Донецкий национальный технический университет Украина, 83001, Донецк, ул. Артема, 66. Применение технологии локальной экстраполяции для высокоточного решения задачи Коши на SIMD-структурах Abstract Feldman L.P., Nazarova I.A. Application of the local extrapolation technology for high exact solution of Cauchy’s problem on the SIMD-structures. Parallel algorithms for system’s of ordinary differential equations on base the method of local extrapolation provided high exact solutions are considerated. The potential system and algorithm parallelism is exploited. Obtained algorithms are realized on SIMD-structures with matrix topology. The estimations of the execution time, acceleration and efficiency parallel solution are defined. Введение В различных предметных областях часто встречаются в качестве математических моделей изучаемых процессов системы обыкновенных дифференциальных уравнений(СОДУ) большой размерности. Бурное развитие большинства современных естественно-научных направлений поставило проблему создания новых, высокоточных (10 -15 – 10 -20 ) методов численного интегрирования таких систем. Традиционный подход к решению этой проблемы, а именно применение явных методов Рунге-Кутты высоких порядков является неэффективным из-за многократного вычисления правых частей дифференциальных уравнений[1,2]. Методы численного интегрирования, использующие «старшие» производные, также вызывают трудности, поскольку в реальных задачах не всегда можно вычислить производные аналитически. Интересным и плодотворным является подход, базирующийся на идее локальной экстраполяции Ричардсона, привлекательность которого резко возросла в связи с развитием высокопроизводительных параллельных вычислительных систем. 1.Общая характеристика технологии локальной экстраполяции Пусть численно решается задача Коши для системы обыкновенных дифференциальных уравнений с известными начальными условиями: . (1) Интегрирование такой системы с некоторой наперед заданной точностью возможно при существовании механизма управления шагом интегрирования на основе определения апостериорной локальной погрешности. Известными способами решения этой проблемы являются: дублирование шага по правилу Рунге, вложенные методы Рунге-Кутты[1,2]. Метод локальной экстраполяции Ричардсона является обобщением технологии удвоения шага по правилу Рунге. Идея этого метода заключается в многократном измельчении шага интегрирования, и также в многократном применении процесса вычисления, названного локальной экстраполяцией (см. рис. 1) . Рисунок 1 – Технология локальной экстраполяции Решение задачи Коши рассматривается при переходе из точки в точку , — базовая длина шага, . Выбирается ряд натуральных чисел такой, что: и, соответственно, последовательность шагов: , где . Задается опорный численный метод порядка и, выполняя шагов интегрирования длиной , вычисляют приближенное решение исходной задачи: . (2) Выполнив вычисления для ряда последовательных значений i, по реккурентному соотношению(2) определяют экстраполированные значения для произвольных i,j по формуле (3)[]. Этот процесс получил название локальная полиномиальная экстраполяция: . (3) Здесь величина b равна единице в общем случае, в тоже время для симметричных опорных методов, имеющих разложение погрешности по степеням , b равно двум (каждая экстраполяция исключает две степени h вместо одной)2. Достоинство этого метода состоит в том, что он дает полную таблицу результатов вычислений, которые образуют последовательность вложенных методов и позволяют оценить локальную погрешность, выбрать стратегию для методов переменного порядка. Таблица 1 – Экстраполяционная таблица p p+b p+2b p … ……… . … … В таблице (1) — есть приближенное решение задачи Коши, полученное численным методом порядка с шагом . Величина соответствует аппроксимации наивысшего порядка, равного 2k, в случае, если вычислены первые k строк экстраполяционной таблицы, а величина соответствует аппроксимации порядка 2k-2. Для управления шагом интегрирования естественно использовать выражение: . Экстраполяционные методы имеют то преимущество, что у них на каждом шаге можно менять не только длину шага, но и порядок метода. Следовательно, экстраполяцию можно рассматривать как метод с переменным порядком (столбцы таблицы) и переменным шагом интегрирования (строки таблицы). Для получения множества шагов интегрирования используются числовые последовательности, образованные гармоническим рядом, степенями двойки, различными четными рядами чисел. Таким образом, экстраполяционная технология Ричардсона включает численный метод решения задачи Коши, последовательность сеток, рекуррентное правило вычисления значений приближенного решения. Эффективность применения технологии локальной экстраполяции напрямую зависит от правильного выбора и сочетания всех трех составляющих этого метода. 2. Исследование эффективности последовательных алгоритмов на основе метода локальной экстраполяции Наиболее эффективным последовательным методом, реализующим технологию локальной полиномиальной экстраполяции cчитается алгоритм Грэгга-Булирша-Штера(ГБШ), базирующийся на модифицированном методе средней точки [3]: (4) Оценим вычислительную сложность этого метода, используя тот факт, что опорный метод ГБШ имеет второй порядок точности и допускает разложение по степеням . Накладные расходы для получения решения экстраполяционным методом порядка на основе (4) определяются количеством обращений к функции f – правой части уравнения (1): . Несколько числовых рядов положительных чисел применяются для обеспечения последовательности сеток интегрирования[1: 1) = 2,4,8,16,32,64,128,… – ; 2) = 2,4,6,8,10,12,14,16,20,22,… – ; (5) 3) = 2,4,6,8,12,16,24,32,48,64,… – . Эффективность использования той или иной, обязательно четной, последовательности натуральных чисел для расчета шагов интегрирования в методе ГБШ можно оценить по графикам рисунка 2: Рис.2 – Эффективность использования четных последовательностей На рисунке 2 представлены зависимости объемов вычислительной работы от порядка экстраполяционного метода. Наиболее оптимальной последовательностью является вторая четная последовательность с минимальным числом вычислений f, для которой: . Рассмотрим применение рассмотренной выше технологии для решения систем линейных однородных уравнений с постоянными коэффициентами : (6) Для сравнения возьмем оптимальный метод ГБШ и оригинальный экстраполяционный экспоненциальный алгоритм на основе гармонического ряда. Точным решением задачи Коши вида (6) является матричная експонента: , . Приближенное решение (6) можно построить, аппроксимировав матричную экспоненту отрезком ряда Тейлора : , (7) причем, (7) определяет численный метод решения (6) порядка [1]. Время последовательной реализации метода ГБШ без учета локальной экстраполяции: где – время умножения матрицы на вектор; – время умножения вектора на скаляр; – время выполнения сложения векторов. Время последовательной реализации экстраполяционного экспоненциального метода без локальной экстраполяции: , где – время вычисления скалярного произведения матриц; — время сложения двух матриц, – время умножения матрицы на скаляр. Заметим, что для получения метода порядка p в случае ГБШ и любого другого симметричного метода необходимо вычислить строк экстраполяционной таблицы, а для экспоненциального метода – р-1- строку. График зависимости времени обоих последовательных алгоритмов с учетом локальной экстраполяции от порядка метода p и размерности исходной системы линейных дифференциальных уравнений m приведен на рисунках 3 и 4(здесь и далее пунктирная линия соответствует экспоненциальному методу, а сплошная – методу ГБШ). Рисунок 3а — Графики времени выполнения последовательных алгоритмов ГБШ и экспоненциального от размерности системы Рисунок 3б — Графики времени выполнения последовательных алгоритмов ГБШ и экспоненциального от порядка метода Таким образом, экстраполяционный экспоненциальный метод решения СЛОДУ с постоянными коэффициентами даже без наличия — экстраполяции требует существенно меньшего времени для последовательной реализации. 3. Оценка качества параллельных алгоритмов с локальной экстраполяцией Рассмотрим отображение описанных алгоритмических схем с локальной экстраполяцией на многопроцессорные вычислительные системы SIMD-структуры с распределенной памятью. Конфигурацию системы считаем фиксированной: число процессорных элементов и схема их соединения не изменяются в процессе счета. Каждый процессор может выполнить любую арифметическую операцию за один такт, временные затраты, связанные с обращением к памяти отсутствуют. Параллельная реализация рассмотренных методов в применении к СЛОДУ требует распараллеливания таких базовых операций, как матричное и векторное умножение и сложение, поэтому в качестве топологического решения принимается замкнутая решетка или тор. На такой топологической схеме наиболее эффективно выполняются матричные операции. Для простоты изложения рассмотрим случай, когда количество процессорных элементов в строке или столбце матрицы совпадает с размерностью задачи, m. Вычисление матричного умножения и умножения матрицы на вектор может быть выполнено по систолическому алгоритму, который является наиболее эффективным для SIMD- систем[5]: , (8) где , , — времена выполнения одиночных операций умножения, сложения и сдвига. Вычисление систолического умножения матрицы на вектор на базе алгоритма сдваивания, подробно описано в [6] и требует следующего времени выполнения: (9) и . Рисунок 4 – Зависимость времени выполнения параллельных алгоритмов от размерности системы В качестве показателей эффективности параллельных алгоритмов будем использовать коэффициенты ускорения и эффективности. Коэффициент ускорения, получаемого при использовании параллельного алгоритма для p процессоров, по сравнению с последовательным вариантом выполнения вычислений определяется: , как отношение времени решения на однопроцессорной ЭВМ: к времени выполнения параллельного алгоритма: . Величина — параметр, количество входных данных задачи, в данном конкретном случае порядок СОДУ. Величина может быть определена, как минимально возможное время выполнения параллельного алгоритма при использовании неограниченного числа процессоров (концепция паракомпьютера): . Эффективность использования параллельным алгоритмом процессоров при решении задачи определяется : . Величина эффективности определяет среднюю долю времени выполнения алгоритма, в течении которой процессоры реально используются для решения задачи. Графики зависимости характеристик параллелизма от размерности исходной задачи приведены на рисунках 6 и 7. При расчете характеристик эффективности считаем, что последовательные варианты этих схем реализованы на однопроцессорной ВС с быстродействием арифметического процессора и объемом ОЗУ равным суммарному объему всех ЗУ арифметических процессоров и с необходимым числом внешних устройств, имеющих скорости обмена такие же, что и МПВС типа SIMD. Рис. 6 – Зависимость коэффициентов ускорения от порядка системы Рис. 7 – Зависимость коэффициентов эффективности от порядка системы Таким образом, предложенный экспоненциальный метод со схемой локальной экстраполяции имеет высокие характеристики потенциального и реального параллелизма при решении СЛОДУ. Определение характеристик параллелизма осуществлялось с помощью пакета MathematicaÒ (Wolfram Research Inc.). Заключение Численный эксперимент и проведенный сравнительный анализ различных последовательных и параллельных вычислительных схем на основе технологии локальной экстраполяции показал как достоинства этого подхода, так и его недостатки. Основным недостатком этого подхода является его высокая вычислительная сложность при решении задач общего вида, даже для наиболее эффективных из известных последовательных методов, таких как ГБШ. Областью приложения экстраполяционных методов остаются высокоточные применения, где такие накладные расходы могут быть оправданными. Одним из путей решения этой проблемы является разработка новых узкоспециализированных методов в рамках общей технологии экстраполяции, а также построение параллельных алгоритмов и схем их отображения на структуры многопроцессорных ВС. Перспективным направлением дальнейших исследований является оценка области и характера устойчивости полученных методов, анализ влияния различных топологий соединения процессоров на характеристики качества параллельных алгоритмов, разработка параллельных алгоритмов для решения СОДУ общего вида. источники: http://meshok.net/item/238826587_%D0%9E_%D0%90%D1%80%D1%83%D1%88%D0%B0%D0%BD%D1%8F%D0%BD_%D0%A7%D0%B8%D1%81%D0%BB%D0%B5%D0%BD%D0%BD%D1%8B%D0%B5_%D1%80%D0%B5%D1%88%D0%B5%D0%BD%D0%B8%D1%8F_%D0%BE%D0%B1%D1%8B%D0%BA%D0%BD%D0%BE%D0%B2%D0%B5%D0%BD%D0%BD%D1%8B%D1%85_%D0%B4%D0%B8%D1%84%D1%84%D0%B5%D1%80%D0%B5%D0%BD%D1%86%D0%B8%D0%B0%D0%BB%D1%8C%D0%BD%D1%8B%D1%85_%D1%83%D1%80%D0%B0%D0%B2%D0%BD%D0%B5%D0%BD%D0%B8%D0%B9_%D0%BD%D0%B0_%D0%A4%D0%BE%D1%80%D1%82%D1%80%D0%B0%D0%BD%D0%B5 http://masters.donntu.org/2009/fvti/scheglov/library/feldman3/index.htm Вам также может понравиться Решение систем уравнений методом постановки Решение задач по математике онлайн //mailru,yandex,google,vkontakte,odnoklassniki,instagram,wargaming,facebook,twitter,liveid,steam,soundcloud,lastfm Уравнения двух пересекающихся плоскостей онлайн Линия пересечения плоскостей онлайн С помощю этого онлайн калькулятора можно найти линию пересечения Второе слагаемое в уравнении бернулли Лекция 4 4.1. Уравнение Бернулли для жидкости Рассмотрим поток жидкости, проходящий по трубопроводу переменно­го Вычислить корень уравнения с дробями Решить уравнение с дробями онлайн При помощи калькулятора можно решать уравнение с дробями.