| Автор cайта: Владимир Потемкин fortran-90@yandex.ru |
 |
Обыкновенные дифференциальные уравнения
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
Y=(/1.0, 0.0, 0.0, 0.5/)
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
Вернуться к оглавлению Скачать DE10
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
Y=(/1.0, 0.0, 0.0, 0.5/)
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
Вернуться к оглавлению Скачать DE15
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.
Вернуться к оглавлению Скачать DE19
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,
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.
j=0
Полученное скорректированное значение принимается за значение функции в данном узле и программа переходит к вычислениям в следующем узле и так продолжается, пока не будет достигнута конечная точка интегрирования.
Для упрощения вычислений неявная формула Адамса может быть преобразована к виду
k-1
yn+1(c) = yn+1(p) - αk-1•h•∑ ∇jfn + αk-1•h•fn+1.
j=0
Коэффициенты αk и βk для явной и неявной формул Адамса в разностной форме для разного порядка аппроксимации k можно найти в книге [Б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
Y=(/1.0, 0.0, 0.0, 0.5/)
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
Вернуться к оглавлению Скачать DE20
DE21
Программа DE21 решает систему из m обыкновенных дифференциальных уравнений первого порядка yi'(x) = fi(x, y1(x), y2(x), ..., ym(x)), i=1, 2, ..., m методом Адамса четвертого порядка аппроксимации с автоматическим выбором величины шага интегрирования. Предполагается, что решение системы имеет непрерывные производные до пятого порядка включительно.
В каждом узле интегрирования решение отыскивается по схеме прогноз-коррекция, которая изложена ранее в описании к программе DE20.
В каждом узле интегрирования программа вычисляет погрешность полученного решения. Если погрешность оказалась меньше заданного значения ε, программа переходит к следующему узлу. Если погрешность больше ε, программа уменьшает шаг в два раза и повторяет вычисления в данном узле с уменьшенной величиной шага. В том случае, когда погрешность решения оказывается меньше, чем ε/32 несколько шагов подряд, шаг интегрирования удваивается. Вычисления продолжаются, пока не будет достигнута конечная точка промежутка.
Вызов программы
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) - массив, который на входе должен содержать начальные значения функций в точке A, на выходе будет содержать решение системы в точке 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
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)
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
Вернуться к оглавлению Скачать DE21
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
Y(1)=1.0; DY(1)=0.0
Y(2)=0.0; DY(2)=0.5
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
Вернуться к оглавлению Скачать DE30
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].
Вызов программы
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.
Вернуться к оглавлению Скачать DE34
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
Y(1)=1.0; DY(1)=0.0
Y(2)=0.0; DY(2)=0.5
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
Вернуться к оглавлению Скачать DE35
