DE10
Программа DE10 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Рунге-Кутта 4-го порядка аппроксимации с постоянным шагом без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 5 порядка включительно.
Решение системы вычисляется по формулам:
k1 = h•f(x0, y0)
k2 = h•f(x0 + h/2, y0 + k1/2)
k3 = h•f(x0 + h/2, y0 + k2/2)
k4 = h•f(x0 + h, y0 + k3)
y1 = y0 + (k1 + 2k2 + 2k3 + k4)/6
Вызов программы
call DE10(DS, A, B, N, Y)
Параметры программы
DS, A, B, N - входные параметры;
Y - входной и выходной параметр;
DS(X, Y, DY) - процедура, которая вычисляет значения производных в точке X:
subroutine DS(x, Y, DY) real, intent(in):: X, Y(m) real:: DY(m) DY(1) = <функция № 1 от X, Y(1), Y(2), ..., Y(m)> . . . . . . . . . . . . . . . . . . . . DY(m) = <функция № m от X, Y(1), Y(2), ..., Y(m)> return end subroutine DS
Real A, B - начальная и конечная точки интегрирования;
Integer N - число промежутков интегрирования от точки A до точки B. N должно быть >=1;
Real Y(1:m) - массив, который на входе должен содержать начальные значения функций в точке A, на выходе будет содержать решение системы в точке B.
Пример
! Интегрирование системы дифференциальных уравнений
! первого порядка методом Рунге-Кутта.
program TestDE10
use NML
implicit none
integer, parameter:: m=4
real:: A, B, X, Y(m), E(m)
integer:: i, n, cnt=0
!begin
A=0.0; B=4.0; n=256
! Начальные значения в точка A
Y=(/1.0, 0.0, 0.0, 0.5/)
! Вычисление точных значений в точке B
E(1)=exp(-B)+B
E(2)=-exp(-B)+1.0
E(3)=0.5*B*exp(2.0*B)
E(4)=0.5*exp(2.0*B)+B*exp(2.0*B)
! Решение системы методом Рунге-Кутта
call DE10(DS, A, B, n, Y)
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
20 format(/(' E',I1,' =',F14.6,' Y',I1,' =',F14.6))
30 format(/' cnt =',I5)
contains
subroutine DS(X, Y, DY)
real, intent(in):: X, Y(:)
real:: DY(:)
!begin
DY(1)=Y(2)
DY(2)=Y(2)+2.0*Y(1)-4.0*Y(3)*exp(-2.0*X)-1.0
DY(3)=Y(4)
DY(4)=2.0*Y(4)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1 ! Счётчик вычисления функции
return
end subroutine DS
end program TestDE10
E1 = 4.018316 Y1 = 4.018196
E2 = 0.981684 Y2 = 0.981648
E3 = 5961.916016 Y3 = 5960.825684
E4 = 13414.310547 Y4 = 13407.820312
cnt = 1024
DE15
Программа DE15 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Рунге-Кутта 5-го порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:
k1 = h•f(x0, y0)
k2 = h•f(x0 + h/4, y0 + k1/4)
k3 = h•f(x0 + 3h/8, y0 + 3k1/32 + 9k2/32)
k4 = h•f(x0 + 12h/13, y0 + 1932k1/2197 - 7200k2/2197 + 7296k3/2197)
k5 = h•f(x0 + h, y0 + 439k1/316 - 8k2 + 3680k3/513 - 845k4/4104)
k6 = h•f(x0 + h/2, y0 - 8k1/27 + 2k2 - 3544k3/2565 + 1859k4/4104 - 11k5/40)
y1 = y0 + 16k1/135 + 6656k3/12825 + 28561k4/56430 - 9k5/50 + 2k6/55
На каждом шаге интегрирования системы программа вычисляет асимптотическую оценку погрешности решения:
e = k1/360 - 128k3/4275 - 2197k4/75240 + k5/50 + 2k6/55
Если погрешность оказывается меньше заданного значения ε, программа переходит к следующему шагу, если больше - программа повторяет вычисления с уменьшенной величиной шага. Величина нового шага hnew расчитывается по формуле:
5 ______
hnew = hold•√ ε/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| меньше минимально-допустимой величины шага hmin, которая вычисляется программой. Параметр H приобретает значение hmin;
Error=2, если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины εmin, которая вычисляется программой. Параметр Eps приобретает значение εmin;
Error=65, если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps. Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
Пример
! Интегрирование системы дифференциальных уравнений
! первого порядка методом Рунге-Кутта
! с автоматическим выбором величины шага.
program TestDE15
use NML
implicit none
integer, parameter:: m=4
real:: A, B, H, Eps, Y(m), E(m)
integer:: i, Error, cnt=0
!begin
A=0.0; B=4.0
Eps=1.0E-7
H=0.03125
! Начальные значения в точка A
Y=(/1.0, 0.0, 0.0, 0.5/)
! Вычисление точных значений в точке B
E(1)=exp(-B)+B
E(2)=-exp(-B)+1.0
E(3)=0.5*B*exp(2.0*B)
E(4)=0.5*exp(2.0*B)+B*exp(2.0*B)
! Решение системы методом Рунге-Кутта
call DE15(DS, A, B, H, Eps, Y, Error)
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt, H
20 format(/(' E',I1,' =',F14.6,' Y',I1,' =',F14.6))
30 format(/' cnt =',I5,' Hout =',F9.6)
contains
subroutine DS(X, Y, DY)
real, intent(in):: X, Y(:)
real:: DY(:)
!begin
DY(1)=Y(2)
DY(2)=Y(2)+2.0*Y(1)-4.0*Y(3)*exp(-2.0*X)-1.0
DY(3)=Y(4)
DY(4)=2.0*Y(4)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine DS
end program TestDE15
E1 = 4.018316 Y1 = 4.018269
E2 = 0.981684 Y2 = 0.981577
E3 = 5961.916016 Y3 = 5961.760254
E4 = 13414.310547 Y4 = 13412.780273
cnt = 606 Hout = 0.039385
DE19
Программа DE19 вычисляет значение производной f0 и разности назад ∇f0, ∇2f0, ..., ∇kf0 порядка k в точке x0 (т.н. фронт метода), необходимые для начала интегрирования системы дифференциальных уравнений первого порядка по методу Адамса. Для нахождения этих разностей применяется алгоритм разгона, который основан на итерационном применении явных формул Адамса с последовательно увеличивающимся порядком аппроксимации. Используя начальное значение y0 в точке x0 и формулу Адамса первого порядка, вычисляем
f0 = f(x0, y0), y1(1) = y0 + hf0, f1(1) = f(x1, y1(1)).
Находим ∇f1(1) = f1(1) - f0 и полагаем ∇f0 = ∇f1(1). Используя формулу Адамса второго порядка, вычисляем
y1(2) = y0 + h(f0 + (1/2)∇f0), y2(2) = y1(2) + h(f1(2) + (1/2)∇f1(2)).
Находим f2(2), ∇f2(2), ∇2f2(2) и полагаем
∇2f0 = ∇2f1(2) = ∇2f2(2),
∇f0 = ∇2f2(2) - 2∇2f2(2).
Используя формулу Адамса третьего порядка, вычисляем
y1(3) = y0 + h(f0 + (1/2)∇f0 + (5/12)∇2f0),
y2(3) = y1(3) + h(f1(3) + (1/2)∇f1(3) + (5/12)∇2f1(3)),
y3(3) = y2(3) + h(f2(3) + (1/2)∇f2(3) + (5/12)∇2f2(3)).
Находим f3(3), ∇f3(3), ∇2f3(3), ∇2f3(3) и полагаем
∇3f0 = ∇3f1 = ∇3f2 = ∇3f3(3),
∇2f0 = ∇3f3(3) - 3∇3f3(3),
∇f0 = ∇f3(3) - 3∇2f3(3) + 3∇3f3(3).
Последовательно применяя формулы Адамса все более высокого порядка получаем все нужные нам разности ∇f0, ∇2f0, ..., ∇kf0 до заданного порядка 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
Программа DE20 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Адамса с порядком аппроксимации от 1 до 6 без контроля точности. Предполагается, что решение системы имеет непрерывные производные на один порядок больше применяемого порядка аппроксимации.
В каждом узле интегрирования решение отыскивается по схеме прогноз-коррекция, которая заключается в следующем. Исходя из точки xn, в которой задано значение функции yn и фронт метода fn, ∇fn, ∇2fn, ..., ∇k-1fn, вычисляется решение в точке xn+1=xn+h по явной формуле Адамса (прогноз):
k-1
yn+1(p) = yn + h•∑ αj∇jfn, (1)
j=0
где ∇j обозначает разность назад j-го порядка, k - порядок аппроксимации. По предсказанному значению yn+1(p) в точке xn+1 вычисляется значение производной fn+1 = f(xn+1, yn+1(p)) и новый фронт метода ∇fn+1, ∇2fn+1, ..., ∇k-1fn+1. Теперь предсказанное значение yn+1(p) можно уточнить по неявной формуле Адамса (коррекция):
k-1
yn+1(c) = yn + h•∑ βj∇jfn+1. (2)
j=0
Полученное скорректированное значение принимается за значение функции в данном узле и программа переходит к вычислениям в следующем узле и так продолжается, пока не будет достигнута конечная точка интегрирования.
Итерационный процесс (1),(2) может быть преобразован к виду [Б6]:
k-1
yn+1(c) = yn+1(p) - αk-1•h•∑ ∇jfn + αk-1•h•fn+1.
j=0
Коэффициенты αk и βk для явной и неявной формул Адамса в разностной форме для разного порядка аппроксимации k приведены в таблице 1 ниже. Увеличение поядка аппроксимации осуществляется простым добавлением новых слогаемых в формулах (1) и (2). Локальная погрешность равна коэффициенту при первом отброшенном члене суммы. В таблице 2 приведены формулы Адамса в классической форме.
| j | α | β |
|---|---|---|
| 0 | 1 | 1 |
| 1 | 1/2 | - 1/2 |
| 2 | 5/12 | - 1/12 |
| 3 | 3/8 | - 1/24 |
| 4 | 251/720 | - 19/720 |
| 5 | 95/288 | - 3/160 |
| 6 | 19087/60480 | - 863/60480 |
| Пор | Формулы | Локаль.погр |
|---|---|---|
| Явные формулы | ||
| 1 | yn+1 = yn + hfn | (1/2)h2y(2) |
| 2 | yn+2 = yn+1 + (1/2)h(3fn+1 - fn) | (5/12)h3y(3) |
| 3 | yn+3 = yn+2 + (1/12)h(23fn+2 - 16fn+1 + 5fn) | (3/8)h4y(4) |
| 4 | yn+4 = yn+3 + (1/24)h(55fn+3 - 59 fn+2 + 37fn+1 - 9fn) | (251/720)h5y(5) |
| 5 | yn+5 = yn+4 + (1/720)h(1901fn+4 - 2774fn+3 + 2616fn+2 - 1274fn+1 - 251fn) | (95/288)h6y(6) |
| Неявные формулы | ||
| 1 | yn+1 = yn + hfn+1 | -(1/2)h2y(2) |
| 1 | yn+1 = yn + (1/2)h(fn+1 + fn) | -(1/12)h3y(3) |
| 2 | yn+2 = yn+1 + (1/12)h(5fn+2 + 8fn+1 = fn) | -(1/24)h4y(4) |
| 3 | yn+3 = yn+2 + (1/24)h(9fn+3 + 19fn+2 - 5fn+1 + fn) | -(19/720)h5y(5) |
| 4 | yn+4 = yn+3 + (1/720)h(251fn+4 + 646fn+3 - 264fn+2 + 106fn+1 - 19fn) | -(3/160)h6y(6) |
Вызов программы
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.
Пример
! Решение системы дифференциальных уравнений
! первого порядка методом Адамса
program TestDE20
use NML
implicit none
integer, parameter:: m=4, N=256, p=4
real:: A, B, Y(m), E(m)
integer:: i, Error, cnt=0
!begin
A=0.0; B=4.0
! Начальные значения в точка A
Y=(/1.0, 0.0, 0.0, 0.5/)
! Вычисление точных значений в точке B
E(1)=exp(-B)+B
E(2)=-exp(-B)+1.0
E(3)=0.5*B*exp(2.0*B)
E(4)=0.5*exp(2.0*B)+B*exp(2.0*B)
! Решение системы методом Адамса
call DE20(DS, A, B, N, m, p, Y)
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
20 format(/(' E',I1,' =',F14.6,' Y',I1,' =',F14.6))
30 format(/' cnt =',I5)
contains
subroutine DS(X, Y, DY)
real, intent(in):: X, Y(:)
real:: DY(:)
!begin
DY(1)=Y(2)
DY(2)=2.0*Y(1)+Y(2)-4.0*Y(3)*exp(-2.0*X)-1.0
DY(3)=Y(4)
DY(4)=2.0*Y(4)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine DS
end program TestDE20
E1 = 4.018316 Y1 = 4.018270
E2 = 0.981684 Y2 = 0.981151
E3 = 5961.916016 Y3 = 5962.983398
E4 = 13414.310547 Y4 = 13417.371094
cnt = 267
DE21
Программа DE21 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Адамса четвертого порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до пятого порядка включительно.
Решение системы осуществляется по схеме прогноз-коррекция, которая изложена ранее в описании к программе DE20.
В каждом узле интегрирования программа вычисляет погрешность полученного решения. Если погрешность оказалась меньше заданного значения ε, программа переходит к следующему узлу. Если погрешность больше ε, программа уменьшает шаг в два раза и повторяет вычисления в данном узле с уменьшенной величиной шага. В том случае, когда погрешность решения оказывается меньше, чем ε/32 несколько шагов подряд, шаг интегрирования удваивается. Вычисления продолжаются, пока не будет достигнута конечная точка промежутка.
Погрешность вычисляется на основе анализа асимптотических формул остаточных членов явного и неявного методов Адамса.
ρ(p) = y(e) - y(p) = (251/720)·h5·y(v)(ξ) ≈ (251/720)·h·∇4fn+1
ρ(c) = y(e) - y(c) = -(19/720)·h5·y(v)(ξ) ≈ -(19/720)·h·∇4fn+1
где ρ(p) и ρ(c) - погрешности явного и неявного методов, y(e) - точное решение (неизвестное), ξ - точка из промежутка (xn, xn+1). Вычитая из первой формулы вторую, имеем:
y(c) - y(p) = (251/720 + 19/720)·h5·y(v)(ξ) = (270/720)·h5·y(v)(ξ),
следовательно h5·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·∇4fn+1 = -(19/720)·h·∇4fn+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| меньше минимально-допустимой величины шага hmin, которая вычисляется программой. Параметр H приобретает значение hmin;
Error=2, если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины εmin, которая вычисляется программой. Параметр Eps приобретает значение εmin;
Error=3, если интегрирование нельзя начать, т.к. погрешность решения, определенная программой при вычислении начальных значений, меньше Eps. Чтобы начать интегрирование, нужно или увеличить величину шага H, или уменьшить значение Eps;
Error=65, если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps. Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
В процессе вычислений программа DE21 вызывает программу DE19.
Пример
!* Решение системы дифференциальных уравнений первого порядка
!* методом Адамса с автоматическим выбором величины шага интегрирования
program TestDE21
use NML
implicit none
integer, parameter:: m=4
real:: A, B, H, Eps, Y(m), E(m)
integer:: i, Error, cnt=0
!begin
A=0.0; B=-4.0
H=0.001953125 !H=1/2**9
Eps=0.5E-7
! Начальные значения в точка A
Y(1)=exp(-A)+A
Y(2)=-exp(-A)+1.0
Y(3)=0.5*A*exp(2.0*A)
Y(4)=0.5*exp(2.0*A)+A*exp(2.0*A)
! Решение системы методом Адамса
call DE21(DS, A, B, H, Eps, Y, Error)
! Вычисление точных значений в точке B
E(1)=exp(-B)+B
E(2)=-exp(-B)+1.0
E(3)=0.5*B*exp(2.0*B)
E(4)=0.5*exp(2.0*B)+B*exp(2.0*B)
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, Error, H, A, cnt
20 format(/(' E',I1,' =',F14.6,' Y',I1,' =',F14.6))
30 format(/' Error =',I2,' H =',E14.6,' X =',E14.6,' cnt ='I5)
contains
subroutine DS(X, Y, DY)
real, intent(in):: X, Y(:)
real:: DY(:)
!begin
DY(1)=Y(2)
DY(2)=Y(2)+2.0*Y(1)-4.0*Y(3)*exp(-2.0*X)-1.0
DY(3)=Y(4)
DY(4)=2.0*Y(4)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine DS
end program TestDE21
E1 = 50.598148 Y1 = 50.598289
E2 = -53.598148 Y2 = -53.598419
E3 = -0.000671 Y3 = -0.000671
E4 = -0.001174 Y4 = -0.001174
Error = 0 H = -0.625000E-01 X = -0.400000E+01 cnt = 118
DE30
Программа DE30 решает систему из m обыкновенных дифференциальных уравнений второго порядка yi"(x) = fi(x, y1(x), ..., ym(x), y1'(x), ..., ym'(x)), i=1, 2, ..., m методом Рунге-Кутта 4-го порядка аппроксимации без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:
k1 = h•f(x0, y0, y0')
k2 = h•f(x0 + ½h, y0 + ½hy0', y0' + ½k1)
k3 = h•f(x0 + ½h, y0 + ½hy0' + ¼hk2, y0' + ½k2)
k4 = h•f(x0 + h, y0 + hy0' + ½hk2, y0' + k3)
y1 = y0 + hy0' + h(k1 + k2 + k3)/6
y1' = y0' +(k1 + 2k2 + 2k3 + k4)/6
Вызов программы
call DE30(DS2, A, B, N, Y, DY)
Параметры программы
DS2, A, B, N - входные параметры;
Y, DY - входной и выходной параметр;
DS2(X, Y, DY, DDY) - процедура, которая вычисляет значения вторых производных в точке X:
subroutine DS2(x, Y, DY, DDY) real, intent(in):: X, Y(m) real:: DY(m), DDY(m) DDY(1) = <функция № 1 от X, Y(1), Y(2), ..., Y(m), DY(1), DY(2), ..., DY(m)> . . . . . . . . . . . . . . . . . . . . . . . DDY(m) = <функция № m от X, Y(1), Y(2), ..., Y(m), DY(1), DY(2), ..., DY(m)> return end subroutine DS2
Real A, B - начальная и конечная точки интегрирования;
Integer N - число промежутков интегрирования от точки A до точки B. N должно быть >=1;
Real Y(1:m) - массив, который на входе должен содержать начальные значения функций в точке A, на выходе будет содержать решение системы в точке B;
Real DY(1:m) - массив, который на входе должен содержать начальные значения производных в точке A, на выходе будет содержать значения производных в точке B.
Пример
! Интегрирование системы дифференциальных уравнений
! второго порядка методом Рунге-Кутта.
program TestDE30
use NML
implicit none
integer, parameter:: m=2
real:: A, B, Y(m), DY(m), E(m)
integer:: i, n, cnt=0
!begin
A=0.0; B=4.0; n=256
! Начальные значения в точка A
Y(1)=1.0; DY(1)=0.0
Y(2)=0.0; DY(2)=0.5
! Вычисление точных значений в точке B
E(1)=B+exp(-B)
E(2)=0.5*B*exp(2.0*B)
! Решение системы методом Рунге-Кутта
call DE30(DS2, A, B, n, Y, DY)
print 10, Y, E, cnt
10 format(/' Y1 =',F10.6,' Y2 =',F14.6 &
/' E1 =',F10.6,' E2 =',F14.6 &
/' cnt =',I5)
contains
subroutine DS2(X, Y, DY, DDY)
real, intent(in):: X, Y(:), DY(:)
real:: DDY(:)
!begin
DDY(1)=DY(1)+2.0*Y(1)-4.0*Y(2)*exp(-2.0*X)-1.0
DDY(2)=2.0*DY(2)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine DS2
end program TestDE30
Y1 = 4.018284 Y2 = 5961.837891
E1 = 4.018316 E2 = 5961.916016
cnt = 1024
DE34
Программа DE34 вычисляет значение производной f0, вспомогательной переменной z0 и разности назад ∇f0, ∇2f0, ..., ∇kf0 порядка k в точке x0, необходимые для начала интегрирования системы дифференциальных уравнений второго порядка по методу Штермера. Указанные величины нходятся с помощью итерационного процесса, в котором применяются явная формула Адамса (для нахождения значения производной y0') и аналог фрмулы Адамса для уравнения второго порядка (для нахождения значения переменной z0):
k-1
y1' = y0' + h•∑ αj∇jfn,
j=0
k-1
z1 = y0 + hy0' + h•∑ μj∇jfn,
j=0
y1 = y0 + hz1.
где αj и μj - коэффициенты этих формул, j=0, 1, ..., k. Применяя эти формулы раз за разом c последовательно увеличивающимся порядком к исходным данным x0, y0, y0', мы получаем необходимые значения для начала интегрирования системы по методу Штермера из точки x0. Более подробную информацию об этой процедуре можно найти в книге [Б6].
Явная и неявная фомулы Штермера:
yn+1 = yn + ∇yn + h2·(fn + (1/12)·∇2fn + (1/12)·∇3fn + (19/240)·∇4fn + (3/40)·∇5fn + ...)
yn+1 = yn + ∇yn + h2·(fn+1 - ∇fn+1 + (1/12)·∇2fn+1 - (1/240)·∇4fn+1 - (1/240)·∇5fn+1 + ...)
Вызов программы
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
Программа DE35 решает систему из m обыкновенных дифференциальных уравнений второго порядка yi"(x) = fi(x, y1(x), ..., ym(x), y1'(x), ..., ym'(x)), i=1, 2, ..., m методом Штермера с порядком аппроксимации от 1 до 6 без контроля точности. Предполагается, что решение системы имеет непрерывные производные на два порядка больше применяемого порядка аппроксимации.
В каждом узле интегрирования решение отыскивается по схеме прогноз-коррекция, которая заключается в следующем. В точке xn+1=xn+h вычисляется приближенное решение по явной формуле Штермера (прогноз):
k-1
yn+1(p) = yn + ∇yn + h2•∑ γj∇jfn,
j=0
где ∇j обозначает разность назад j-го порядка, k - порядок аппроксимации. Полученное приближение можно уточнить, выполнив вычисления по неявной формуле Штермера (коррекция):
k-1
yn+1(с) = yn + ∇yn + h2•∑ δj∇jfn+1.
j=0
Здесь γj и δj обозначают коэффициенты явной и неявной формул Штермера. С целью уменьшения вычислительной погрешности при интегрировании системы вводится новая переменная zn и формулы Штермера пиобретают вид:
k-1
zn+1(p) = zn + h•∑ γj∇jfn,
j=0
k-1
zn+1(c) = zn + h•∑ δj∇jfn+1,
j=0
yn+1(c) = yn + hzn+1(c).
Значения первой производной, также как и в программе DE20, вычисляются по схеме прогноз-коррекция по явной и неявной формулам Адамса:
k-1
yn+1'(p) = yn' + h•∑ αj∇jfn,
j=0
k-1
yn+1'(c) = yn' + h•∑ βj∇jfn+1.
j=0
Для упрощения вычислений формулы коррекции слегка видоизменяются, чтобы в их состав входили уже вычисленные предсказанные значения zn+1(p) и yn+1'(p). Окончательно вычислительный процесс выглядит следующим образом. Сначала по явным формулам вычисляются значения функции и производной в точке xn+1:
k-1
zn+1(p) = zn + h•∑ γj∇jfn,
j=0
yn+1(p) = yn + hzn+1(p),
k-1
yn+1'(p) = yn' + h•∑ αj∇jfn.
j=0
Далее вычисляется значение fn+1=f(xn+1, yn+1(p), yn+1'(p)) и новый фронт метода. Скорректированные значения вычисляются по формулам:
k-1
zn+1(c) = zn+1(p) - γk-1•h•∑ ∇jfn + γk-1•h•fn+1.
j=0
yn+1(c) = yn + hzn+1(c),
k-1
yn+1'(c) = yn+1'(p) - αk-1•h•∑ ∇jfn + αk-1•h•fn+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.
Пример
! Решение системы дифференциальных уравнений
! второго поядка методом Штермера
program TestDE35
use NML
implicit none
integer, parameter:: m=2, N=256, p=4
real:: A, B, Y(m), DY(m), E(m)
integer:: i, cnt=0
!begin
A=0.0; B=4.0
! Начальные значения в точка A
Y(1)=1.0; DY(1)=0.0
Y(2)=0.0; DY(2)=0.5
! Вычисление точных значений в точке B
E(1)=B+exp(-B)
E(2)=0.5*B*exp(2.0*B)
! Решение системы методом Штермера
call DE35(DS2, A, B, N, m, p, Y, DY)
print 10, Y, E, cnt
10 format(/' Y1 =',F10.6,' Y2 =',F14.6 &
/' E1 =',F10.6,' E2 =',F14.6 &
/' cnt =',I5)
contains
subroutine DS2(X, Y, DY, DDY)
real, intent(in):: X, Y(:), DY(:)
real:: DDY(:)
!begin
DDY(1)=DY(1)+2.0*Y(1)-4.0*Y(2)*exp(-2.0*X)-1.0
DDY(2)=2.0*DY(2)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine DS2
end program TestDE35
Y1 = 4.018204 Y2 = 5961.664062
E1 = 4.018316 E2 = 5961.916016
cnt = 267
DE23
Программа DE23 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом трапеций (неявным методом Адамса 1 порядка аппроксимации):
yn+1' = yn + (1/2)h(f(xn+1, yn+1) + f(xn, yn))
Для повышения устойчивости итерационный процесс для нахождения yn+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 итерации не удалось достигнуть заданной точности по какой-либо компоненте решения.
Пример 1
program Trap_1
! решение системы уравнений 1-го порядка
! y' = -20y + z
! z' = 19y - 2z
! методом трапеций с постоянным шагом
use NML
implicit none
integer, parameter:: m=2, N=256, NW=0
character, parameter:: met='T'
real:: Eps
real:: A, B, H, Y(m), E(m)
integer:: NE
integer:: cnt=0, cnt2=0
real Em1, Em21
integer i
!begin
Eps=1.0E-6
A=0.0; B=1.0
! Начальные значения в точка A
Y=(/2.0, 18.0/)
! Вычисление точных значений в точке B
Em1=exp(-1.0); Em21=exp(-21.0)
E(1)=Em1+Em21
E(2)=19.0*Em1-Em21
! Решение системы методом трапеций
call DE23R(F, FJ, A, B, Y, N, Eps, met, NW, NE)
print 19, N, Eps
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
print 31, cnt2
print 32, NE
19 format(/(4X,'N =',I6,' Eps =',F11.7))
20 format(/(4X,'E',I1,' =',E14.7,' Y',I1,' =',E14.7))
30 format(/4X,'Number of function calculations =',I5)
31 format(4X,'Number of calculations of the Jacobian =',I5)
32 format(4X,'Number of error limits exceeded =',I5)
contains
subroutine F(X, Y, DY)
! вычисление производных в точке Х
real, intent(in):: X, Y(:)
real:: DY(:)
DY(1)=-20.0*Y(1)+Y(2)
DY(2)=19.0*Y(1)-2.0*Y(2)
cnt=cnt+1
return
end subroutine F
subroutine FJ(X, Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: X, Y(:)
real:: YJ(:,:)
YJ(1,1)=-20.0; YJ(2,1)=19.0;
YJ(1,2)=1.0; YJ(2,2)=-2.0;
cnt2=cnt2+1
return
end subroutine FJ
end program Trap_1
N = 256 Eps = 0.0000010
E1 = 0.3678795E+00 Y1 = 0.3678789E+00
E2 = 0.6989709E+01 Y2 = 0.6989700E+01
Number of function calculations = 513
Number of calculations of the Jacobian = 1
Number of error limits exceeded = 0
http://num-anal.srcc.msu.ru/lib_na/cat/de/de23r.htm
Пример 2
program Trap_2
! решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! методом трапеций с постоянным шагом
use NML
implicit none
integer, parameter:: m=4, N=2048, NW=1
character, parameter:: met='T'
real:: Eps
real:: A, B, H, Y(m), E(m)
integer:: NE
integer:: cnt=0, cnt2=0
real C0, C1, Lam
integer i
!begin
Eps=1.0E-6
A=0.0; B=4.0
! Начальные значения в точка A
Y=(/1.0, 0.0, 0.0, 0.5/)
! Вычисление точных значений в точке B
E(1)=exp(-B)+B
E(2)=-exp(-B)+1.0
E(3)=0.5*B*exp(2.0*B)
E(4)=0.5*exp(2.0*B)+B*exp(2.0*B)
! Решение системы методом трапеций
call DE23R(F, FJ, A, B, Y, N, Eps, met, NW, NE)
print 19, N, Eps
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
print 31, cnt2
print 32, NE
19 format(/(4X,'N =',I6,' Eps =',F11.7))
20 format(/(4X,'E',I1,' =',E14.7,' Y',I1,' =',E14.7))
30 format(/4X,'Number of function calculations =',I5)
31 format(4X,'Number of calculations of the Jacobian =',I5)
32 format(4X,'Number of error limits exceeded =',I5)
contains
subroutine F(X, Y, DY)
! вычисление производных в точке Х
real, intent(in):: X, Y(:)
real:: DY(:)
DY(1)=Y(2)
DY(2)=2.0*Y(1)+Y(2)-4.0*Y(3)*exp(-2.0*X)-1.0
DY(3)=Y(4)
DY(4)=2.0*Y(4)+(Y(1)-X)*exp(3.0*X)
cnt=cnt+1
return
end subroutine F
subroutine FJ(X, Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: X, Y(:)
real:: YJ(:,:)
YJ(1,1)=0.0; YJ(2,1)=2.0; YJ(3,1)=0.0; YJ(4,1)=exp(3.0*X)
YJ(1,2)=1.0; YJ(2,2)=1.0; YJ(3,2)=0.0; YJ(4,2)=0.0
YJ(1,3)=0.0; YJ(2,3)=-4.0*exp(-2.0*X); YJ(3,3)=0.0; YJ(4,3)=0.0
YJ(1,4)=0.0; YJ(2,4)=0.0; YJ(3,4)=1.0; YJ(4,4)=2.0
cnt2=cnt2+1
end subroutine FJ
end program Trap_2
N = 2048 Eps = 0.0000010
E1 = 0.4018316E+01 Y1 = 0.4019418E+01
E2 = 0.9816844E+00 Y2 = 0.9850047E+00
E3 = 0.5961916E+04 Y3 = 0.5963318E+04
E4 = 0.1341431E+05 Y4 = 0.1344188E+05
Number of function calculations = 4097
Number of calculations of the Jacobian = 2048
Number of error limits exceeded = 0
DE26
Программа DE23 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m по формуле дифференцирования назад 2-го поядка аппроксимации:
yn+2' = (4/3)yn+1 - (1/3)yn + (2/3)hf(xn+2, yn+2)
Для повышения устойчивости итерационный процесс для нахождения yn+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 итерации не удалось достигнуть заданной точности по какой-либо компоненте решения.
Пример 1
program Difback
! решение уравнения
! y'=-100*(y-sin(x))
! y(0)=1
!
use NML
implicit none
integer, parameter:: m=1, N=256, NW=0
real:: Eps
real:: A, B, H, Y(m), E(m)
integer:: NE
integer:: cnt=0, cnt2=0
real C0, C1, Lam
integer i
!begin
Eps=1.0E-6
A=0.0; B=3.14159265358979
! Начальные значения в точка A
Y=(/1.0/)
! Вычисление точных значений в точке B
Lam=100.0; C1=Lam/(1+Lam*Lam); C0=1.0+C1
E(1)=C0*exp(-Lam*B)-C1*(-Lam*sin(B)+cos(B))
! Решение системы по формуле диф.назад
call DE26R(F, FJ, A, B, Y, N, Eps, NW, NE)
print 19, N, Eps
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
print 31, cnt2
print 32, NE
19 format(/(4X,'N =',I6,' Eps =',F11.7))
20 format(/(4X,'E',I1,' =',E14.7,' Y',I1,' =',E14.7))
30 format(/4X,'Number of function calculations =',I5)
31 format(4X,'Number of calculations of the Jacobian =',I5)
32 format(4X,'Number of error limits exceeded =',I5)
contains
subroutine F(X, Y, DY)
! вычисление производных в точке Х
real, intent(in):: X, Y(:)
real:: DY(:)
DY(1)=-Lam*(Y(1)-sin(X))
cnt=cnt+1
return
end subroutine F
subroutine FJ(X, Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: X, Y(:)
real:: YJ(:,:)
YJ(1,1)=-100.0
cnt2=cnt2+1
end subroutine FJ
end program Difback
N = 256 Eps = 0.0000010
E1 = 0.9998913E-02 Y1 = 0.9999443E-02
Number of function calculations = 491
Number of calculations of the Jacobian = 2
Number of error limits exceeded = 0
Пример 2
program Difback_2
! решение системы уравнений 1-го порядка
! y' = -20y + z
! z' = 19y - 2z
!
use NML
implicit none
integer, parameter:: m=2, N=256, NW=0
real:: Eps
real:: A, B, H, Y(m), E(m)
integer:: NE
integer:: cnt=0, cnt2=0
real Em1, Em21
integer i
!begin
Eps=1.0E-6
A=0.0; B=1.0
! Начальные значения в точка A
Y=(/2.0, 18.0/)
! Вычисление точных значений в точке B
Em1=exp(-1.0); Em21=exp(-21.0)
E(1)=Em1+Em21
E(2)=19.0*Em1-Em21
! Решение системы по формуле диф.назад
call DE26R(F, FJ, A, B, Y, N, Eps, NW, NE)
print 19, N, Eps
print 20, (i, E(i), i, Y(i), i=1,m)
print 30, cnt
print 31, cnt2
print 32, NE
19 format(/(4X,'N =',I6,' Eps =',F11.7))
20 format(/(4X,'E',I1,' =',E14.7,' Y',I1,' =',E14.7))
30 format(/4X,'Number of function calculations =',I5)
31 format(4X,'Number of calculations of the Jacobian =',I5)
32 format(4X,'Number of error limits exceeded =',I5)
contains
subroutine F(X, Y, DY)
! вычисление производных в точке Х
real, intent(in):: X, Y(:)
real:: DY(:)
DY(1)=-20.0*Y(1)+Y(2)
DY(2)=19.0*Y(1)-2.0*Y(2)
cnt=cnt+1
return
end subroutine F
subroutine FJ(X, Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: X, Y(:)
real:: YJ(:,:)
YJ(1,1)=-20.0; YJ(2,1)=19.0;
YJ(1,2)=1.0; YJ(2,2)=-2.0;
cnt2=cnt2+1
end subroutine FJ
end program Difback_2
N = 256 Eps = 0.0000010
E1 = 0.3678795E+00 Y1 = 0.3678854E+00
E2 = 0.6989709E+01 Y2 = 0.6989823E+01
Number of function calculations = 347
Number of calculations of the Jacobian = 2
Number of error limits exceeded = 0
DE53
Программа DE53 решает задчу Коши для жесткой автономной системы из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом типа Розенброка 3-го порядка аппроксимации с автоматическим выбором величины шага интегрирования. В случае неавтономной системы y' = f(x, y(x)) введением дополнительной переменной y'm+1 = 1, ym+1(x0) = x0 она приводится к автономному виду.
Решение системы вычисляется по формулам:
Dnk1 = h•f(yn)
Dnk2 = h•f(yn + β21k1)
Dnk3 = h•f(yn + β31k1 + β32k2)
yn+1 = yn + p1k1 + p2k2 + p3k3
где Dn = E + α•h•fn′;
E — единичная матрица; h — шаг интегрирования;
fn′ = ∂f(yn)/∂y — матрица Якоби системы;
α, β21, β31, β32, p1, p2, p3 — числовые коэффициенты, определяющие свойства точности и устойчивости алгоритма.
В каждом узле интегрирования программа вычисляет погрешность полученного решения. Для оценки погрешности используется идея вложенных методов, а именно, дополнительно вычисляется
yn+1/2 = yn + b1k1 + b2k2
Теперь оценить погрешность можно по фомуле
εn(2) ≈ (1/c)·Dn-1·(yn+1 + yn+1/2)
где c - коэффициент, Dn-1 - матрица, обратная к матрице Dn.
В программе вычисляется коэффициент q2 из формулы
q23·║εn(2)║ = C·ε
где ε - заданная пользователем точность вычислений, C - коэффициент. Если q2 < 1, то программа повторяет вычисления из прежней точки с меньшей величиной шага. Если q2 > 1, то программа считает, что требуемая точность получена и переходит к вычислению следующего шага. В обоих случаях новый шаг расчитывается по формуле hnew = k·q2·hold, где k - понижающий коэффициент.
Норма ║ξ║ расчитывается по формуле
║ξ║ = max1≤j≤m{|ξj|/(|ynj| + r)}
где ynj - j-ая компонента решения, r – малый положительный параметр. Если |ynj| < r, применяется абсолютная погрешность, иначе относительная погрешность.
[Б17]
Вызов программы
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| меньше минимально-допустимой величины шага hmin, которая вычисляется программой. Параметр H приобретает значение hmin;
Error=2, если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины εmin, которая вычисляется программой. Параметр Eps приобретает значение εmin;
Error=65, если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps. Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
В процессе вычислений программа DE53 вызывает программу AIG2R.
Пример 1
program Rosenbrock1
! решение автономной системы уравнений 1-го порядка
! y' = -20y + z
! z' = 19y - 2z
! методом Розенброка 3-го порядка
use NML
implicit none
integer, parameter:: m=2, nW=0
real:: Eps
real:: X, Xout, H, H0, E(m)
real:: Y(m)
integer:: nE
integer:: cnt=0, cnt2=0
real Em1, Em21
integer i
!begin
Eps=1.0E-4 ! Погрешность
H=1.0/7168.0; H0=H ! Начальный шаг
X=0.0; Xout=1.0 ! Начальная и конечная точки
! Начальные значения в точке Х
Y=(/2.0, 18.0/)
! Вычисление точных значений в точке Xout
Em1=exp(-1.0*Xout); Em21=exp(-21.0*Xout)
E(1)=Em1+Em21
E(2)=19.0*Em1-Em21
! Решение системы методом Розенброка
call DE53R(F, FJ, X, Xout, H, Y, Eps, nW, nE)
print 20, (i, E(i), i, Y(i), i=1,m)
print 29, H0, H
print 30, cnt
print 31, cnt2
print 32, nE
20 format(/(4X,'E',I1,' =',F14.6,' Y',I1,' =',F14.6))
29 format(/4X,'The initial step = ',E15.7/4X,'The step at the exit of the program = ',E15.7)
30 format(4X,'Number of calculations of the function = ',I5)
31 format(4X,'Number of calculations of the Jacobian = ',I5)
32 format(4X,'Error = ',I2)
stop
contains
subroutine F(Y, DY)
! вычисление производных в точке Х
real, intent(in):: Y(:)
real:: DY(:)
DY(1)=-20.0*Y(1)+Y(2)
DY(2)=19.0*Y(1)-2.0*Y(2)
cnt=cnt+1
return
end subroutine F
subroutine FJ(Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: Y(:)
real:: YJ(:,:)
YJ(1,1)=-20.0; YJ(2,1)=19.0;
YJ(1,2)=1.0; YJ(2,2)=-2.0;
cnt2=cnt2+1
end subroutine FJ
end program Rosenbrock1
E1 = 0.367879 Y1 = 0.367849
E2 = 6.989709 Y2 = 6.989136
The initial step = 0.1395089E-03
The step at the exit of the program = 0.1391107E-03
Number of calculations of the function = 21567
Number of calculations of the Jacobian = 1
Error = 0
Пример 2
program Rosenbrock2
! решение автономной системы уравнений 1-го порядка
! y1' = -2.0*y1*y1*(y2*y2 + y3*y3)
! y2' = y2 + y3
! y3' = -y2 + y3
! методом Розенброка 3-го порядка
use NML
implicit none
integer, parameter:: m=3, nW=1
real:: Eps
real:: X, Xout, H, H0, E(m)
real:: Y(m)
integer:: nE
integer:: cnt=0, cnt2=0
integer i
!begin
Eps=1.0E-4 ! Погрешность
H=1.0/8192.0; H0=H ! Начальный шаг
X=0.0; Xout=3.0 ! Начальная и конечная точки
! Начальные значения в точке Х
Y=(/1.0, 0.0, 1.0/)
! Вычисление точных значений в точке Xout
E(1)=exp(-2.*Xout)
E(2)=exp(Xout)*sin(Xout)
E(3)=exp(Xout)*cos(Xout)
! Решение системы методом Розенброка
call DE53R(F, FJ, X, Xout, H, Y, Eps, nW, nE)
print 20, (i, E(i), i, Y(i), i=1,m)
print 29, H0, H
print 30, cnt
print 31, cnt2
print 32, nE
20 format(/(4X,'E',I1,' =',F14.6,' Y',I1,' =',F14.6))
29 format(/4X,'The initial step = ',E15.7/4X,'The step at the exit of the program = ',E15.7)
30 format(4X,'Number of calculations of the function = ',I5)
31 format(4X,'Number of calculations of the Jacobian = ',I5)
32 format(4X,'Error = ',I2)
stop
contains
subroutine F(Y, DY)
! вычисление производных в точке Х
real, intent(in):: Y(:)
real:: DY(:)
DY(1)=-2.*Y(1)**2*(Y(2)**2+Y(3)**2)
DY(2)=Y(2)+Y(3)
DY(3)=-Y(2)+Y(3)
cnt=cnt+1
return
end subroutine F
subroutine FJ(Y, YJ)
! вычисление матрицы Якоби в точке Х
real, intent(in):: Y(:)
real:: YJ(:,:)
YJ(1,1)=-4.*Y(1)*(Y(2)**2+Y(3)**2); YJ(2,1)=0.; YJ(3,1)=0.
YJ(1,2)=-4.*Y(1)**2*Y(2); YJ(2,2)=1.; YJ(3,2)=-1.
YJ(1,3)=-4.*Y(1)**2*Y(3); YJ(2,3)=1.; YJ(3,3)=1.
cnt2=cnt2+1
return
end subroutine FJ
end program Rosenbrock2
E1 = 0.2478752E-02 Y1 = 0.2476111E-02
E2 = 0.2834471E+01 Y2 = 0.2842367E+01
E3 = -0.1988453E+02 Y3 = -0.1989177E+02
The initial step = 0.1220703E-03
The step at the exit of the program = 0.1390711E-03
Number of calculations of the function = 64752
Number of calculations of the Jacobian = 21582
Error = 0
Пример 3
program Rose3
! решение не автономного уравнения
! y'=-100*(y-sin(x))
! y(0)=1
! методом Розенброка 3-го порядка
use NML
implicit none
integer, parameter:: m=2, nW=1
real:: Eps
real:: X, Xout, H, H0, E(m)
real:: Y(m)
integer:: nE
integer:: cnt=0, cnt2=0
real C0, C1, Lam
integer i
!begin
Eps=1.0E-4 ! Погрешность
H=1.0/8192.0; H0=H ! Начальный шаг
X=0.0; Xout=3.14159265358979 ! Начальная и конечная точки
! Начальные значения в точке Х
Y=(/1.0, 0.0/)
! Вычисление точных значений в точке Xout
Lam=100.0; C1=Lam/(1+Lam*Lam); C0=1.0+C1
E(1)=C0*exp(-Lam*Xout)-C1*(-Lam*sin(Xout)+cos(Xout))
E(2)=Xout
! Решение системы методом Розенброка
call DE53R(F, FJ, X, Xout, H, Y, Eps, nW, nE)
print 20, (i, E(i), i, Y(i), i=1,m)
print 29, H0, H
print 30, cnt
print 31, cnt2
print 32, nE
20 format(/(4X,'E',I1,' =',F14.6,' Y',I1,' =',F14.6))
29 format(/4X,'The initial step = ',E15.7/4X,'The step at the exit of the program = ',E15.7)
30 format(4X,'Number of calculations of the function = ',I5)
31 format(4X,'Number of calculations of the Jacobian = ',I5)
32 format(4X,'Error = ',I2)
stop
contains
subroutine F(Y, DY)
! вычисление производных
real, intent(in):: Y(:)
real:: DY(:)
DY(1)=-Lam*(Y(1)-sin(X))
DY(2)=1.0
cnt=cnt+1
return
end subroutine F
subroutine FJ(Y, YJ)
! вычисление матрицы Якоби
real, intent(in):: Y(:)
real:: YJ(:,:)
YJ(1,1)=-100.0; YJ(2,1)=0.0;
YJ(1,2)=cos(Y(2)); YJ(2,2)=0.0;
cnt2=cnt2+1
end subroutine FJ
end program Rose3
E1 = 0.009999 Y1 = 0.009948
E2 = 3.141593 Y2 = 3.142450
The initial step = 0.1220703E-03
The step at the exit of the program = 0.1376044E-03
Number of calculations of the function = 68181
Number of calculations of the Jacobian = 22726
Error = 0
DT10
Программа DT10 решает краевую задачу для линейного дифференциального уравнения 2-го порядка y"(x) + q(x)y'(x) + p(x)y(x) = f(x) на промежутке [a, b] с граничными условиями y(a)=ya, y(b)=yb методом прогонки. Первая и вторая производные заменяются разностными схемами y' = (yi+1 - yi-1)/2h и y" = (yi+1 - 2yi + yi-1)/h², i = 1, ..., n-1; x0 = a, xn = b; h = xi+1 - xi - шаг сетки. Подстановка прводит к системе из n-2 уравнений:
ai·yi-1 + bi·yi + ci·yi+1 = di, i = 1, ..., n-1
где ai = 2 - h·q(xi), bi = -2(2 - h²·p(xi)), ci = 2 + h·q(xi), di = 2h²·f(xi). Два недостающих уравнения беруться из граничных условий y(x0) = y0 и y(xn) = yn.
Получившаяся система из 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 задает разрядность сетки.
Пример
program test_dt10
! Краевая задача для уравнения
! y" - 2x·y' - 2·y = 4x
! y(0) = 1; y(1) = ℮ - 1
use NML
implicit none
integer, parameter:: n=21
real X(n), Y(n), E(n)
real A, B, YA, YB
real h, R
integer i
!begin
A=0.0; B=1.0
! начальные значения
YA=exp(A*A)-A; YB=exp(B*B)-B
! точные значения
h=(B-A)/real(n-1)
do i=1, n
R=A+h*real(i-1); E(i)=exp(R*R)-R
end do
! X(1)=A; X(n)=B
! решение задачи методом прогонки
call DT10R(Q, P, F, A, B, YA, YB, X, Y)
print 10
print 11, (X(i), Y(i), E(i), i=1,n,2)
10 format(/10X,'X',14X,'Y',17X,'E'/)
11 format(4X,F10.6,4X,E14.7,4X,E14.7)
stop
contains
real function Q(X)
real, intent(in):: X
Q=-2.0*X
end function Q
real function P(X)
real, intent(in):: X
P=-2.0
end function P
real function F(X)
real, intent(in):: X
F=4.0*X
end function F
end program test_dt10
X Y E
0.000000 0.1000000E+01 0.1000000E+01
0.100000 0.9101487E+00 0.9100502E+00
0.200000 0.8409864E+00 0.8408108E+00
0.300000 0.7944074E+00 0.7941743E+00
0.400000 0.7737821E+00 0.7735109E+00
0.500000 0.7843144E+00 0.7840254E+00
0.600000 0.8336133E+00 0.8333294E+00
0.700000 0.9325696E+00 0.9323162E+00
0.800000 0.1096676E+01 0.1096481E+01
0.900000 0.1348016E+01 0.1347908E+01
1.000000 0.1718282E+01 0.1718282E+01
DT11
Программа DT11 решает краевую задачу для линейного дифференциального уравнения 2-го порядка y"(x) + q(x)y'(x) + p(x)y(x) = f(x) на промежутке [a, b] со смешанными граничными условиями методом прогонки. Первая и вторая производные заменяются разностными схемами y' = (yi+1 - yi-1)/2h и y" = (yi+1 - 2yi + yi-1)/h², i = 1, ..., n-1; x0 = a, xn = b; h = xi+1 - xi - шаг сетки. Подстановка прводит к системе из n-2 уравнений:
ai·yi-1 + bi·yi + ci·yi+1 = di, i = 1, ..., n-1
где ai = 2 - h·q(xi), bi = -2(2 - h²·p(xi)), ci = 2 + h·q(xi), di = 2h²·f(xi). Два недостающих уравнения берём из смешанных граничных условий:
y'(x0) + α·y(x0) = α1
y'(xn) + β·y(xn) = β1
Данные граничные условия заменяем разностными схемами:
(1/2h)·(y1 - y-1) + α·y0 = α1
(1/2h)·(yn+1 - yn-1) + β·yn = β1
из которых выражаем неизвестные y-1 и yn+1:
y-1 = y1 - 2h·(α1 - α·y0)
yn+1 = yn-1 + 2h·(β1 - β·yn)
Подставляем их в уравнения:
a0·y-1 + b0·y0 + c0·y1 = d0
an·yn-1 + bn·yn + cn·yn+1 = dn
и получаем два недостающих уравнения для решения системы:
(b0 + 2h·α·a0)·y0 + (a0 + c0)·y1 = d0 + 2h·α1·a0
(an + cn)·yn-1 + (bn - 2h·β·cn)·yn = dn - 2h·β1·cn
Получившаяся система из 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
program test_dt11
! Краевая задача для уравнения
! y" + 2x/(x² + 1)·y' + 2/(x² + 1)·y = 4x/(x² + 1)³
! y'(-1) - y(-1) = 0
! y'(1) + y(1) = 0
use NML
implicit none
integer, parameter:: n=81
real X(n), Y(n), E(n)
real A, B, Alpha, Alpha1, Beta, Beta1
real h, R
integer i
!begin
A=-1.0; B=1.0
! начальные значения
Alpha=-1.0; Alpha1=0.0
Beta=1.0; Beta1=0.0
! точные значения
h=(B-A)/real(n-1)
do i=1, n
R=A+h*real(i-1); E(i)=1.0/(R*R+1)
end do
! X(1)=A; X(n)=B
! решение задачи методом прогонки
call DT11R(Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1, X, Y)
print 10
print 11, (X(i), Y(i), E(i), i=1,n,8)
10 format(/10X,'X',14X,'Y',17X,'E'/)
11 format(4X,F10.6,4X,E14.7,4X,E14.7)
stop
contains
real function Q(X)
real, intent(in):: X
Q=2.0*X/(X*X+1)
end function Q
real function P(X)
real, intent(in):: X
P=2.0/(X*X+1)
end function P
real function F(X)
real, intent(in):: X
F=4.0*X*X/(X*X+1)**3
end function F
end program test_dt11
X Y E
-1.000000 0.4998799E+00 0.5000000E+00
-0.800000 0.6096135E+00 0.6097562E+00
-0.600000 0.7351372E+00 0.7352942E+00
-0.400000 0.8619106E+00 0.8620690E+00
-0.200000 0.9613886E+00 0.9615385E+00
0.000000 0.9998574E+00 0.1000000E+01
0.200000 0.9613931E+00 0.9615385E+00
0.400000 0.8619165E+00 0.8620689E+00
0.600000 0.7351429E+00 0.7352940E+00
0.800000 0.6096182E+00 0.6097561E+00
1.000000 0.4998826E+00 0.5000000E+00
Пример 2
program boundary_task_2
! Краевая задача для уравнения
! y" - 2x·y' - 2·y = 4x
! y'(0) + 3·y(0) = 2
! y'(1) + (4℮ - 3)/(℮ - 1)·y(1) = 6℮ - 4
use NML
implicit none
integer, parameter:: n=20
real X(0:n), Y(0:n), E(0:n)
! размерность массивов X и Y 21 элемент
real A, B, Alpha, Alpha1, Beta, Beta1
real h, R, e1
integer i
!begin
A=0.0; B=1.0
! начальные значения
Alpha=3.0; Alpha1=2.0
e1=exp(1.0)
Beta=(4.0*e1-3.0)/(e1-1.0); Beta1=6.0*e1-4.0
! точные значения
h=(B-A)/real(n)
do i=0, n
R=A+h*real(i); E(i)=exp(R*R)-R
end do
! X(0)=A; X(n)=B
! решение задачи методом прогонки
call DT11R(Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1, X, Y)
print 10
print 11, (X(i), Y(i), E(i), i=0,n,2)
10 format(/10X,'X',14X,'Y',17X,'E'/)
11 format(4X,F10.6,4X,E14.7,4X,E14.7)
stop
contains
real function Q(X)
real, intent(in):: X
Q=-2.0*X
end function Q
real function P(X)
real, intent(in):: X
P=-2.0
end function P
real function F(X)
real, intent(in):: X
F=4.0*X
end function F
end program boundary_task_2
X Y E
0.000000 0.1000236E+01 0.1000000E+01
0.100000 0.9102049E+00 0.9100502E+00
0.200000 0.8408605E+00 0.8408108E+00
0.300000 0.7940893E+00 0.7941743E+00
0.400000 0.7732548E+00 0.7735109E+00
0.500000 0.7835505E+00 0.7840254E+00
0.600000 0.8325718E+00 0.8333294E+00
0.700000 0.9311924E+00 0.9323162E+00
0.800000 0.1094882E+01 0.1096481E+01
0.900000 0.1345695E+01 0.1347908E+01
1.000000 0.1715279E+01 0.1718282E+01