Автор cайта: Владимир Потемкин fortran-90@yandex.ru
Обыкновенные дифференциальные уравнения
DE10
Программа DE10 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Рунге-Кутта 4-го порядка аппроксимации с постоянным шагом без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 5 порядка включительно.
Решение системы вычисляется по формулам:
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 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Рунге-Кутта 5-го порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:
Если погрешность оказывается меньше заданного значения ε, программа переходит к следующему шагу, если больше - программа повторяет вычисления с уменьшенной величиной шага. Величина нового шага 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 вычисляет значение производной f0 и разности назад ∇f0, ∇2f0, ..., ∇kf0 порядка k в точке x0 (т.н. фронт метода), необходимые для начала интегрирования системы дифференциальных уравнений первого порядка по методу Адамса. Для нахождения этих разностей применяется алгоритм разгона, который основан на итерационном применении явных формул Адамса с последовательно увеличивающимся порядком аппроксимации. Используя начальное значение y0 в точке x0 и формулу Адамса первого порядка, вычисляем
Последовательно применяя формулы Адамса все более высокого порядка получаем все нужные нам разности ∇f0, ∇2f0, ..., ∇kf0 до заданного порядка k включительно. Более подробную информацию об алгоритме разгона можно найти в книге [Б6].
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 обыкновенных дифференциальных уравнений первого порядка 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 приведены формулы Адамса в классической форме.
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 решает систему из 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 решает систему из m обыкновенных дифференциальных уравнений второго порядка yi"(x) = fi(x, y1(x), ..., ym(x), y1'(x), ..., ym'(x)), i=1, 2, ..., m методом Рунге-Кутта 4-го порядка аппроксимации без контроля точности. Предполагается, что решение системы имеет непрерывные производные до 6 порядка включительно.
Решение системы вычисляется по формулам:
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 вычисляет значение производной 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].
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 обыкновенных дифференциальных уравнений второго порядка 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 решает систему из 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 DE23R(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
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 решает систему из 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 можно оставить обращение матрицы только перед первой итерацией, а в остальных итерациях пользоваться уже вычисленным значением. Это позволяет значительно уменьшить объём вычислений. В некоторых случаях (для систем уравнений с постоянными коэффициентами) матрица Якоби не зависит от аргумента и функций, т.е. является числовой матрицей. Тогда матрица Якоби обращается только один раз в начале выполнения программы и далее на всём выполнении программы используется уже вычисленное значение.
В процессе вычислений программа DE26 вызывает программу AIG2R обращения матрицы.
Вызов программы
call DE26R(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 решает задчу Коши для жесткой автономной системы из 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 она приводится к автономному виду.
где 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
Пакет программ для решения системы обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Адамса/Гира с заданным начальным условием yi(x0) = yi0, i=1, 2, ..., m. DDRIV2 содержит опции, позволяющие решать как жесткие, так и нежесткие дифференциальные уравнения. Программа производит вычисления с переменными двойной точности.
T - Независимая переменная. На входе при первом обращении значение T равно начальной точке интегрирования, на выходе - в точке T находится вычисленное решение.
Y - Вектор зависимых переменных. На входе при первом обращении в Y задаются начальные значения. На выходе Y содержит вектор вычисленного решения. Массив Y передаётся через список параметров в подпрограмму F и подпрограмму-функцию G, которые должен написать пользователь. Информация, необходимая для F и G, может быть размещена в этом массиве начиная с (N+1)-й компоненты. (Замечание. Изменения, вносимые пользователем в первые Y компонент, появятся лишь в случае повторного обращения к программе, т.е. после изменения MSTATE на +1(-1).)
F - Подпрограмма, которую должен написать пользователь. Её имя должно быть описано в операторе EXTERNAL вызывающей программы пользователя. Это подпрограмма следующего вида:
SUBROUTINE F(N, T, Y, YDOT)
DOUBLE PRECISION Y(*), YDOT(*)
...
YDOT(1)=...
...
YDOT(N)=...
END SUBROUTINE F
Она вычисляет YDOT = F(Y,T) - правые части дифференциальных уравнений. Здесь Y - массив длины не меньше N. Фактический размер массива Y описывается пользователем в программе, вызывающей DDRIV2. При обращении к подпрограмме F первые N элементов массива Y содержат вычисляемые приближения к компонентам решения. Пользователь не должен изменять значения элементов в этом массиве. В данной подпрограмме YDOT - массив дины N. Пользователю следует вычислить YDOT(I) только для I от 1 до N. Как правило, возврат из F снова передаёт управление подпрограмме DDRIV2. Однако, если пользователь желает прервать вычисления, т.е. передать управление в программу, вызвавшую DDRIV2, ему следует присвоить N нулевое значение. DDRIV2 сигнализирует об этом, возвратив MSTATE со значением 6(-6). Смена значения N в F не повлияет на значение параметра N, указанного в списке параметров DDRIV2.
TOUT - На входе: Точка, в которой требуется вычислить решение.
MSTATE - Указатель целого типа, описывающий режим интегрирования. Пользователь может задать MSTATE равным +1 или -1. Если MSTATE положительно, то подпрограмма, миновав при интегрировании точку TOUT, интерполирует решение. Это наиболее эффективный режим. Если MSTATE отрицательно, то подпрограмма будет подбирать шаги интегрирования чтобы попасть точно в точку TOUT. Этот режим потребуется в случае сингулярности правее точки TOUT.
Смысл модулей значений MSTATE таков:
1. (На входе) Первое обращение к подпрограмме. Данное значение должно быть задано пользователем. При всех последующих обращениях пользователю необходимо лишь следить за величиной MSTATE. До тех пор, пока к DDRIV2 не обратились заново, пользователем может быть изменён только знак MSTATE. (Для удобства пользователя, который хотел бы отказаться от задания начальных условий, DDRIV2 допускает вызов со значениями параметров MSTATE = +1 (-1) и TOUT=T. В таком случае подпрограмма вернёт MSTATE т.е. MSTATE = +1 (-1).)
2. (На выходе) Успешное завершение интегрирования. Если требуется продолжить интегрирование (в том же направлении), то просто увеличьте TOUT и обратитесь снова. Все прочие параметры будут установлены автоматически.
3. (На выходе) Неудачное завершение. Выполнено 1000 шагов интегрирования, но точка TOUT. не достигнута. Пользователь может продолжить интегрирование, просто снова обратившись к DDRIV2. Если это не ошибка в постановке задачи, то наиболее вероятной причиной может быть попытка решить жёсткую систему уравнений в режиме нежёсткого интегрирования. (См. ниже описание параметра MINT.)
4. (На выходе) Неудачное завершение. Затребована слишком высокая точность. EPS было увеличено до уровня, рассматриваемого программой как разумный. Пользователь может продолжить интегрирование, просто снова обратившись к DDRIV2.
5. (На выходе) Корень был найден в точке, лежащей левее точки TOUT (правее, если T > TOUT). Пользователь может продолжить интегрирование к точке TOUT, просто снова обратившись к DDRIV2.
6. (На выходе) Неудачное завершение. В подпрограмме F для N было установлено значение 0.
7. (На выходе) Неудачное завершение. В подпрограмме-функции G для N было установлено значение 0. (См. ниже описание G.)
NROOT - Число уравнений, для которых ищется корень. Если NROOT есть нуль, то режим поиска корня не активируется. Этот режим полезен, когда требуется получить результат в точках, неизвестных заранее, но зависящих от решения, например в точках, где некоторая компонента решения принимает определённое значение. Поиск корня ведётся с помощью подпрограммы-функции G, написанной пользователем (см. ниже описание G). DDRIV2 пытается найти значение T, при котором одно из уравнений меняет знак. DDRIV2 может найти в лучшем случае один корень в уравнении на внутреннем шаге интегрирования, а затем возвратит решение или в точке TOUT или в корне, в зависимости от того, какая из этих точек встретится первой в направлении интегрирования. Номер уравнения, корень которого был найден, записывается в элемент массива IWORK.
Замечание. NROOT никогда не изменяется в DDRIV2.
EPS - На входе: требуемая относительная точность во всех компонентах решения. Значение EPS = 0 допускается. На выходе: откорректированная относительная точность если входное значение оказалось слишком малым. Значение EPS следует задать - по возможности - в пределах разумного - большим, поскольку объём вычислений в DDRIV2 возрастает при уменьшении EPS.
EWT - (На входе) Нуль задачи, т.е. наименьшее значение (компоненты) решения, которое может быть допущено по смыслу задачи. Эта величина используется внутри DDRIV2 при вычислении массива YWT(I) = MAX(ABS(Y(I),EWT)). Оценка погрешности на шаге, делённая на YWT(I), поддерживается меньшей, чем EPS. При EWT равным нулю имеем просто-напросто контроль относительной погрешности. Однако, задав EWT слишком малым можно столкнуться с увеличением времени счёта.
MINT - (На входе) Указатель метода интегрирования.
MINT = 1 - используются методы Адамса; рекомендуется для нежёстких задач.
MINT = 2 - используются «жёсткие» методы Гира (т.е. формулы дифференцирования назад); рекомендуется для жёстких задач.
MINT = 3 - программа динамически выбирает или методы Адамса, если задача оказывается нежёсткой, или методы Гира, когда она является жёсткой.
Значение MINT не может быть переопределено без дополнительного обращения к программе, т.е. без задания для переменной MSTATE значения 1.
MITER - (На входе) Индикатор итерационного метода.
MITER = 0 - Означает функциональную итерацию. Это значение рекомендуется для нежёстких систем уравнений.
MITER = 1 - Означает метод хорд с аналитическим якобианом. В этом случае пользователь должен предоставить подпрограмму JACOBN (см. описание ниже).
MITER = 2 - Означает метод хорд с якобианом, вычисляемым внутренне с помощью конечных разностей.
MITER = 3 - Означает метод хорд с поправками, вычисляемыми пользователями с помощью подпрограммы USERS (см. описание USERS ниже). Эта опция позволяет пользователю принимать все решения по матричной алгебре и хранению данных. При использовании значения MITER = 3 выполнение подпрограммы FA не требуется, даже если значение IMPL не равно 0. Дополнительную информацию об использовании этой опции смотрите в разделе IV-E ниже.
MITER = 4 - Означает то же самое, что и MITER = 1, но предполагается, что матрицы A и Якоби являются ленточными.
MITER = 5 - Означает то же самое, что и MITER = 2, но предполагается, что матрицы A и Якоби являются ленточными.
Значение MITER не должно изменяться во время выполнения программы.
IMPL - (На входе) Индикатор неявного метода.
IMPL = 0 - Означает решение задачи dY(I)/dT = F(Y(I),T).
IMPL = 1 - Означает решение задачи A*dY(I)/dT = F(Y(I),T), не особенная матрица A (см. описание FA ниже). Для этого параметра допустимы только значения MINT = 1 или 2 и MITER = 1, 2, 3, 4 или 5.
IMPL = 2 - Означает решение определенных систем гибридных дифференциальных/алгебраических уравнений (см. описание FA ниже). Для этого параметра допустимы только значения MINT = 2 и MITER = 1, 2, 3, 4 или 5.
Значение IMPL не должно изменяться во время выполнения программы.
WORK, LENW - (На входе) WORK - вещественный массив длины LENW, используемый внутри программы в качестве рабочего. Пользователь должен отвести в вызывающей программе под этот массив память, например, посредством описания вида:
DOUBLE PRECISION WORK(...)
Длина массива WORK должна быть не меньше, чем
16 × N + 2 × NROOT + 204, если MINT = 1;
N × N + 10 × N + 2 × NROOT + 204, если MINT = 2;
N × N + 17 × N + 2 × NROOT + 204, если MINT = 3;
(установленный размер памяти следует присвоить целой переменной LENW). Содержимое массива WORK не должно изменяться между обращениями к DDRIV2.
IWORK, LENIW - (На входе) IWORK - целый массив длины LENIW, используемый внутри программы в качестве рабочего. Пользователь должен отвести в вызывающей программе под этот массив память, например, посредством описания вида:
INTEGER IWORK(...)
Длина массива WORK должна быть не меньше, чем
21, если MINT = 1;
N + 21, если MINT = 2 или 3;
(установленный размер памяти следует присвоить целой переменной LEIW). Содержимое массива IWORK не должно изменяться между обращениями к DDRIV2.
G - Вещественная подпрограмма-функция, которую должен написать пользователь в том случае, если NROOT не есть 0. В таком случае её имя должно быть описано в операторе EXTERNAL в вызывающей программе. G последовательно вызывается с различными значениями IROOT для того, чтобы получить значение каждой из функций NROOT, корни которых подлежат вычислению. Имеет следующий вид:
double precision function G(N, T, Y, IROOT)
integer:: N, IROOT
double precision:: T, Y(*)
!begin
select case (IROOT)
case (1)
G=...
. . .
case (N)
G=...
end select
end function G
Здесь Y - массив длины не меньшей N; его первые N элементов содержат решение в точке T. Пользователь не должен изменять эти значения. Фактический размер массива Y определяется по его описанию в программе пользователя, вызывающей DDRIV2.
При правильном завершении программы управление возвращается в DDRIV2. Однако, если пользователь желает прервать вычисления, т.е. передать управление в программу, вызвавшую DDRIV2, ему следует присвоить N нулевое значение. DDRIV2 сигнализирует об этом, возвратив MSTATE со значением, равным 7 (-7). В таком случае номер обработанного в G уравнения будет размещён в шестом элементе массива IWORK. Смена значения N в G не повлияет на значение параметра N, указанного в списке параметров DDRIV2.
JACOBN - Подпрограмма, которую должен написать пользователь когда MITER равно 1 или 4. (См. комментарии к подпрограмме DDRIV3.)
FA - Подпрограмма, которую должен написать пользователь когда IMPL равно 1 или 2 и MITER равно 3. (См. комментарии к подпрограмме DDRIV3.)
USERS - Подпрограмма, которую должен написать пользователь когда MITER равно 3. (См. комментарии к подпрограмме DDRIV3.)
ДОПОЛНИТЕЛЬНАЯ ИНФОРМАЦИЯ ДЛЯ ПОЛЬЗОВАТЕЛЯ
Первые три элемента массива WORK и первые пять элементов массива IWORK будет помещена следующая статистическая информация:
AVGH - среднее значение использованных шагов,
HUSED - последний (успешный) использованный шаг,
AVGORD - среднее значение использованных порядков,
IMXERR - индекс того компонента вектора решения, на котором достигается максимальная ошибка при последней проверке,
NQUSED - последний (успешно) использованный порядок,
NSTEP - число шагов, выполненных после последней инициализации,
NFE - число вычислений правых частей,
NJE -число вычислений матриц-якобианов.
Описание программы SDRIV2 приведено в книге [Б18].
Выражаю искреннюю благодарность Василию Рябову за предоставленные тексты программ.
Пример A
program TestDDRIV2_A
! Решение задачи Коши методом Адамса
use DDRIV
implicit none
! Решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! в промежутке от 0 до 4 с начальными условиями
! y(0)=1, y'(0)=0, z(0)=0, z'(0)=1/2
integer, parameter:: MINT = 1, MITER = 0, NROOT = 0, &
N = 4, LENW = 16*N + 2*NROOT + 204, LENIW = 21
! N is the number of equations
!EXTERNAL F, G, JACOBN, FA, USERS
real(kind=8):: EPS, EWT, T, TOUT, WORK(LENW), Y(N)
integer MSTATE, IWORK(LENIW)
real(kind=8) A, B, H, E(N)
data WORK(1),WORK(2),WORK(3) /3*0.D0/
data IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5) /5*0/
integer i
A = 0.D0 ! Initial point
B = 4.D0 ! Final point, |A|<|B|
H = 4.D0 ! Integration step
Y(1)=exp(-A)+A ! Set initial conditions
Y(2)=-exp(-A)+1.D0
Y(3)=0.5D0*A*exp(2.D0*A)
Y(4)=0.5D0*exp(2.D0*A)+A*exp(2.D0*A)
E(1)=exp(-B)+B ! Calculating the exact values at point B
E(2)=-exp(-B)+1.D0
E(3)=0.5D0*B*exp(2.D0*B)
E(4)=0.5D0*exp(2.D0*B)+B*exp(2.D0*B)
T = A
TOUT = A
MSTATE = 1 ! +1 or -1
EPS = 1.D-5
EWT = 1.D-5
20 call DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT, &
MINT, MITER, WORK, LENW, IWORK, LENIW, G,JACOBN,FA,USERS)
if (MSTATE .GT. 2) go to 30
TOUT = TOUT + H
if (abs(TOUT) .LE. abs(B)) go to 20
30 continue
print 100, MSTATE, TOUT, EPS, EWT
print 101, (i, E(i), i, Y(i), i=1,N)
print 102, WORK(1), WORK(2), WORK(3)
print 103, IWORK(1), IWORK(2), IWORK(3), IWORK(4), IWORK(5)
print 104, IWORK(4)
100 format(/4X,'MSTATE =',I3,' TOUT =',F5.1,' EPS =',D13.5,' EWT =',D13.5)
101 format(/(8X,'E',I1,' =',F15.6,' Y',I1,' =',F15.6))
102 format(/4X,'AVG =',D14.6,'; HUSED =',D14.6,'; AVGORD =',D14.6)
103 format(4X,'IMXERR =',I4,'; NQUSED =',I4,'; NSTEP =',I4,'; NFE =',I4,'; NJE =',I6)
104 format(/8X,'Number of calculations of the function',I5)
contains
subroutine F(N, T, Y, YDOT)
integer:: N
double precision:: T, Y(*)
double precision:: YDOT(*)
!begin
YDOT(1)=Y(2)
YDOT(2)=Y(2)+2.0*Y(1)-4.0*Y(3)*exp(-2.0*T)-1.0
YDOT(3)=Y(4)
YDOT(4)=2.0*Y(4)+(Y(1)-T)*exp(3.0*T)
end subroutine F
double precision function G(N, T, Y, IROOT)
! Необязательная программа.
! Пользователь должен написать её только если NROOT не равно 0.
integer:: N, IROOT
double precision:: T, Y(*)
!begin
select case (IROOT)
case (1)
G=0.D0
case default
G=0.D0
end select
end function G
subroutine JACOBN (N, T, Y, DFDY, MATDIM, ML, MU)
! Необязательная программа вычисления Якобиана.
! Пользователь должен написать её только если MITER равно 1 или 4.
! Если MITER = 1, то частная производня от правой части i-го уравнения
! по j-ой зависимой функции запоминается в элементе DFDY(i,j);
! MATDIM = 1
! Если MITER = 4, то матрица Якоби ленточная и частные производные
! i-го уравнения относительно Y(j) должны быть сохранены в
! элементе DFDY(k,j), где k = i-j+MU+1;
! MATDIM = 2*ML+MU+1; ML и MU - число лент ниже и выше диагонали.
integer:: N, MATDIM, ML, MU
double precision:: T, Y(*)
double precision:: DFDY(MATDIM,*)
!C .
!C .
!C Calculate values of DFDY
!C .
!C .
end subroutine JACOBN
subroutine FA(N, T, Y, A, MATDIM, ML, MU, NDE)
! Необязательная программа.
! Пользователь должен написать её только если
! значение IMPL равно 1 или 2, а значение MITER не равно 3.
! Переменная NDE требуется, если количество уравнений NDE меньше N
integer:: N, MATDIM, ML, MU, NDE
double precision:: T, Y(*)
double precision:: A(MATDIM,*)
!C .
!C .
!C Calculate ALL values of A
!C .
!C .
end subroutine FA
subroutine USERS(Y, YH, YWT, SAVE1, SAVE2, T, H, EL, IMPL, N, NDE, IFLAG)
! Необязательная программа.
! Пользователь должен написать её только если значение MITER равно 3.
integer:: N, NDE, IMPL, IFLAG
double precision Y(*), YH(*), YWT(*), SAVE1(*), SAVE2(*), T, H, EL
!C .
!C
!C
end subroutine USERS
end program TestDDRIV2_A
MSTATE = 2 TOUT = 8.0 EPS = 0.10000D-04 EWT = 0.10000D-04
E1 = 4.018316 Y1 = 4.020487
E2 = 0.981684 Y2 = 0.987607
E3 = 5961.915974 Y3 = 5966.477691
E4 = 13414.310942 Y4 = 13475.233075
AVG = 0.741425D-01; HUSED = 0.149423D+00; AVGORD = 0.448148D+01
IMXERR = 4; NQUSED = 6; NSTEP = 54; NFE = 95; NJE = 0
Number of calculations of the function 95
Пример G
program TestDDRIV2_G
! Решение задачи Коши методом Гира, Якобиан вычисляется в функции JACOBN
use DDRIV
implicit none
! Решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! в промежутке от 0 до 4 с начальными условиями
! y(0)=1, y'(0)=0, z(0)=0, z'(0)=1/2
integer, parameter:: MINT = 2, MITER = 1, NROOT = 0, &
N = 4, LENW = N*N + 10*N + 2*NROOT + 204, LENIW = N + 21
! N is the number of equations
!EXTERNAL F, G, JACOBN, FA, USERS
real(kind=8) EPS, EWT, T, TOUT, WORK(LENW), Y(N)
integer MSTATE, IWORK(LENIW)
real(kind=8) A, B, H, E(N)
data WORK(1),WORK(2),WORK(3) /3*0.D0/
data IWORK(1),IWORK(2),IWORK(3),IWORK(4),IWORK(5) /5*0/
integer i
A = 0.D0 ! Initial point
B = 4.D0 ! Final point, |A|<|B|
H = 4.D0 ! Integration step
Y(1)=exp(-A)+A ! Set initial conditions
Y(2)=-exp(-A)+1.D0
Y(3)=0.5D0*A*exp(2.D0*A)
Y(4)=0.5D0*exp(2.D0*A)+A*exp(2.D0*A)
E(1)=exp(-B)+B ! Calculating the exact values at point B
E(2)=-exp(-B)+1.D0
E(3)=0.5D0*B*exp(2.D0*B)
E(4)=0.5D0*exp(2.D0*B)+B*exp(2.D0*B)
T = A
TOUT = A
MSTATE = 1 ! +1 or -1
EPS = 1.D-5
EWT = 1.D-5
20 call DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT, &
MINT, MITER, WORK, LENW, IWORK, LENIW, G,JACOBN,FA,USERS)
if (MSTATE .GT. 2) go to 30
TOUT = TOUT + H
if (abs(TOUT) .LE. abs(B)) go to 20
30 continue
print 100, MSTATE, TOUT, EPS, EWT
print 101, (i, E(i), i, Y(i), i=1,N)
print 102, WORK(1), WORK(2), WORK(3)
print 103, IWORK(1), IWORK(2), IWORK(3), IWORK(4), IWORK(5)
print 104, IWORK(4)
print 105, IWORK(5)
100 format(/4X,'MSTATE =',I3,' TOUT =',F5.1,' EPS =',D13.5,' EWT =',D13.5)
101 format(/(8X,'E',I1,' =',F15.6,' Y',I1,' =',F15.6))
102 format(/4X,'AVG =',D14.6,'; HUSED =',D14.6,'; AVGORD =',D14.6)
103 format(4X,'IMXERR =',I4,'; NQUSED =',I4,'; NSTEP =',I4,'; NFE =',I4,'; NJE =',I6)
104 format(/8X,'Number of calculations of the function',I5)
105 format(15X,'Number of Jacobian calculations',I5)
contains
subroutine F(N, T, Y, YDOT)
integer:: N
double precision:: T, Y(*)
double precision:: YDOT(*)
!begin
YDOT(1)=Y(2)
YDOT(2)=Y(2)+2.0*Y(1)-4.0*Y(3)*exp(-2.0*T)-1.0
YDOT(3)=Y(4)
YDOT(4)=2.0*Y(4)+(Y(1)-T)*exp(3.0*T)
end subroutine F
double precision function G(N, T, Y, IROOT)
! Необязательная программа.
! Пользователь должен написать её только если NROOT не равно 0.
integer:: N, IROOT
double precision:: T, Y(*)
!begin
select case (IROOT)
case (1)
G=0.D0
case default
G=0.D0
end select
end function G
subroutine JACOBN (N, T, Y, DFDY, MATDIM, ML, MU)
integer:: N, MATDIM, ML, MU
double precision:: T, Y(*)
double precision:: DFDY(MATDIM,*)
DFDY(1,1)=0.D0; DFDY(2,1)=2.D0; DFDY(3,1)=0.D0; DFDY(4,1)=exp(3.D0*T)
DFDY(1,2)=1.D0; DFDY(2,2)=1.D0; DFDY(3,2)=0.D0; DFDY(4,2)=0.D0
DFDY(1,3)=0.D0; DFDY(2,3)=-4.D0*exp(-2.D0*T); DFDY(3,3)=0.D0; DFDY(4,3)=0.D0
DFDY(1,4)=0.D0; DFDY(2,4)=0.D0; DFDY(3,4)=1.D0; DFDY(4,4)=2.D0
end subroutine JACOBN
subroutine FA(N, T, Y, A, MATDIM, ML, MU, NDE)
! Необязательная программа.
! Пользователь должен написать её только если
! значение IMPL равно 1 или 2, а значение MITER не равно 3.
! Переменная NDE требуется, если количество уравнений NDE меньше N
integer:: N, MATDIM, ML, MU, NDE
double precision:: T, Y(*)
double precision:: A(MATDIM,*)
!C .
!C .
!C Calculate ALL values of A
!C .
!C .
end subroutine FA
subroutine USERS(Y, YH, YWT, SAVE1, SAVE2, T, H, EL, IMPL, N, NDE, IFLAG)
! Необязательная программа.
! Пользователь должен написать её только если значение MITER равно 3.
integer:: N, NDE, IMPL, IFLAG
double precision Y(*), YH(*), YWT(*), SAVE1(*), SAVE2(*), T, H, EL
!C .
!C
!C
end subroutine USERS
END PROGRAM TestDDRIV2_G
MSTATE = 2 TOUT = 8.0 EPS = 0.10000D-04 EWT = 0.10000D-04
E1 = 4.018316 Y1 = 4.025814
E2 = 0.981684 Y2 = 1.017153
E3 = 5961.915974 Y3 = 5934.593269
E4 = 13414.310942 Y4 = 13464.759608
AVG = 0.412491D-01; HUSED = 0.542651D-01; AVGORD = 0.376289D+01
IMXERR = 4; NQUSED = 5; NSTEP = 97; NFE = 147; NJE = 15
Number of calculations of the function 147
Number of Jacobian calculations 15
Пример G+
Всё как в предыдущем примере, но Якобиан вычисляется не в функции, а с помощью конечных разностей (MITER = 2). Результаты совпадают, но число вычислений функции увеличилось.
Number of calculations of the function 207
Number of Jacobian calculations 15
Программа 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
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
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
