Обыкновенные дифференциальные уравнения. Библиотека численных методов на языке Fortran 90

DE10

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

k₁ = h•f(x₀, y₀)
k₂ = h•f(x₀ + h/2, y₀ + k₁/2)
k₃ = h•f(x₀ + h/2, y₀ + k₂/2)
k₄ = h•f(x₀ + h, y₀ + k₃)
y₁ = y₀ + (k₁ + 2k₂ + 2k₃ + k₄)/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

Пример

! Интегрирование системы дифференциальных уравнений
! y1'=y2
! y2'=Y2+2*y1-4*y3*exp(-2*x)-1
! y3'=y4
! y4'=2*y4+(y1-x)*exp(3*x)
! с начальными условиями
! y1(0)=1, y2(0)=0, y3(0)=0, y4(0)=0.5
! методом Рунге-Кутта 4-го порядка

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

Вернуться к оглавлению Скачать DE10

DE15

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

k₁ = h•f(x₀, y₀)
k₂ = h•f(x₀ + h/4, y₀ + k₁/4)
k₃ = h•f(x₀ + 3h/8, y₀ + 3k₁/32 + 9k₂/32)
k₄ = h•f(x₀ + 12h/13, y₀ + 1932k₁/2197 - 7200k₂/2197 + 7296k₃/2197)
k₅ = h•f(x₀ + h, y₀ + 439k₁/216 - 8k₂ + 3680k₃/513 - 845k₄/4104)
k₆ = h•f(x₀ + h/2, y₀ - 8k₁/27 + 2k₂ - 3544k₃/2565 + 1859k₄/4104 - 11k₅/40)
y₁ = y₀ + 16k₁/135 + 6656k₃/12825 + 28561k₄/56430 - 9k₅/50 + 2k₆/55

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

e = k₁/360 - 128k₃/4275 - 2197k₄/75240 + k₅/50 + 2k₆/55

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

Вызов программы

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

Параметры программы

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

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

Пример

! Интегрирование системы дифференциальных уравнений
! y1'=y2
! y2'=Y2+2*y1-4*y3*exp(-2*x)-1
! y3'=y4
! y4'=2*y4+(y1-x)*exp(3*x)
! с начальными условиями
! y1(0)=1, y2(0)=0, y3(0)=0, y4(0)=0.5
! методом Рунге-Кутта 5-го порядка
! с автоматическим выбором величины шага

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

Вернуться к оглавлению Скачать DE15

DE19

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

f₀ = f(x₀, y₀), y₁⁽¹⁾ = y₀ + hf₀, f₁⁽¹⁾ = f(x₁, y₁⁽¹⁾).

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

y₁⁽²⁾ = y₀ + h(f₀ + (1/2)∇f₀),
y₂⁽²⁾ = y₁⁽²⁾ + h(f₁⁽²⁾ + (1/2)∇f₁⁽²⁾).

Находим f₂⁽²⁾, ∇f₂⁽²⁾, ∇²f₂⁽²⁾ и полагаем

∇²f₀ = ∇²f₁⁽²⁾ = ∇²f₂⁽²⁾,
∇f₀ = ∇²f₂⁽²⁾ - 2∇²f₂⁽²⁾.

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

y₁⁽³⁾ = y₀ + h(f₀ + (1/2)∇f₀ + (5/12)∇²f₀),
y₂⁽³⁾ = y₁⁽³⁾ + h(f₁⁽³⁾ + (1/2)∇f₁⁽³⁾ + (5/12)∇²f₁⁽³⁾),
y₃⁽³⁾ = y₂⁽³⁾ + h(f₂⁽³⁾ + (1/2)∇f₂⁽³⁾ + (5/12)∇²f₂⁽³⁾).

Находим f₃⁽³⁾, ∇f₃⁽³⁾, ∇²f₃⁽³⁾, ∇²f₃⁽³⁾ и полагаем

∇³f₀ = ∇³f₁ = ∇³f₂ = ∇³f₃⁽³⁾,
∇²f₀ = ∇³f₃⁽³⁾ - 3∇³f₃⁽³⁾,
∇f₀ = ∇f₃⁽³⁾ - 3∇²f₃⁽³⁾ + 3∇³f₃⁽³⁾.

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

Таблица 1
j	α	β
0	1	1
1	1/2	- 1/2
2	5/12	- 1/12
3	3/8	- 1/24
4	251/720	- 19/720
5	95/288	- 3/160
6	19087/60480	- 863/60480

Таблица 2
Пор	Формулы	Локаль.погр
	Явные формулы
1	y_n+1 = y_n + hf_n	(1/2)h²y⁽²⁾
2	y_n+2 = y_n+1 + (1/2)h(3f_n+1 - f_n)	(5/12)h³y⁽³⁾
3	y_n+3 = y_n+2 + (1/12)h(23f_n+2 - 16f_n+1 + 5f_n)	(3/8)h⁴y⁽⁴⁾
4	y_n+4 = y_n+3 + (1/24)h(55f_n+3 - 59 f_n+2 + 37f_n+1 - 9f_n)	(251/720)h⁵y⁽⁵⁾
5	y_n+5 = y_n+4 + (1/720)h(1901f_n+4 - 2774f_n+3 + 2616f_n+2 - 1274f_n+1 - 251f_n)	(95/288)h⁶y⁽⁶⁾
	Неявные формулы
1	y_n+1 = y_n + hf_n+1	-(1/2)h²y⁽²⁾
1	y_n+1 = y_n + (1/2)h(f_n+1 + f_n)	-(1/12)h³y⁽³⁾
2	y_n+2 = y_n+1 + (1/12)h(5f_n+2 + 8f_n+1 - f_n)	-(1/24)h⁴y⁽⁴⁾
3	y_n+3 = y_n+2 + (1/24)h(9f_n+3 + 19f_n+2 - 5f_n+1 + f_n)	-(19/720)h⁵y⁽⁵⁾
4	y_n+4 = y_n+3 + (1/720)h(251f_n+4 + 646f_n+3 - 264f_n+2 + 106f_n+1 - 19f_n)	-(3/160)h⁶y⁽⁶⁾

Вызов программы

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.

Пример

! Решение системы дифференциальных уравнений
! y1'=y2
! y2'=Y2+2*y1-4*y3*exp(-2*x)-1
! y3'=y4
! y4'=2*y4+(y1-x)*exp(3*x)
! с начальными условиями
! y1(0)=1, y2(0)=0, y3(0)=0, y4(0)=0.5
! методом Адамса

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

Вернуться к оглавлению Скачать DE20

DE21

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

Вызов программы

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

Параметры программы

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

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

Пример

! Решение системы дифференциальных уравнений
! y1'=y2
! y2'=Y2+2*y1-4*y3*exp(-2*x)-1
! y3'=y4
! y4'=2*y4+(y1-x)*exp(3*x)
! с начальными условиями
! y1(0)=1, y2(0)=0, y3(0)=0, y4(0)=0.5
! методом Адамса 4-го порядка
! с автоматическим выбором величины шага интегрирования

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

Вернуться к оглавлению Скачать DE21

DE30

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

k₁ = h•f(x₀, y₀, y₀')
k₂ = h•f(x₀ + ½h, y₀ + ½hy₀', y₀' + ½k₁)
k₃ = h•f(x₀ + ½h, y₀ + ½hy₀' + ¼hk₂, y₀' + ½k₂)
k₄ = h•f(x₀ + h, y₀ + hy₀' + ½hk₂, y₀' + k₃)
y₁ = y₀ + hy₀' + h(k₁ + k₂ + k₃)/6
y₁' = y₀' +(k₁ + 2k₂ + 2k₃ + k₄)/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.

Пример

! решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! с начальными условиями
! y(0)=1, y'(0)=0, z(0)=0, z'(0)=0.5
! методом Рунге-Кутта 4-го порядка

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

Вернуться к оглавлению Скачать DE30

DE34

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

Явная и неявная фомулы Штермера:
y_n+1 = y_n + ∇y_n + h²·(f_n + (1/12)·∇²f_n + (1/12)·∇³f_n + (19/240)·∇⁴f_n + (3/40)·∇⁵f_n + ...)
y_n+1 = y_n + ∇y_n + h²·(f_n+1 - ∇f_n+1 + (1/12)·∇²f_n+1 - (1/240)·∇⁴f_n+1 - (1/240)·∇⁵f_n+1 + ...)

Вызов программы

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

Параметры программы

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

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

Вернуться к оглавлению Скачать DE34

DE35

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

Вызов программы

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

Параметры программы

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

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

Пример

! решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! с начальными условиями
! y(0)=1, y'(0)=0, z(0)=0, z'(0)=0.5
! методом Штермера

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

Вернуться к оглавлению Скачать DE35

DE23

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

Вызов программы

call 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
! с начальными условиями y(0)=2, z(0)=18
! методом трапеций с постоянным шагом
use NML
implicit none

    integer, parameter:: m=2, N=256, NW=0
    character:: met
    real:: Eps
    real:: A, B, H, Y(m), E(m)
    integer:: NE
    integer:: cnt=0, cnt2=0
    real Em1, Em21
    integer i
!begin
    met='T'
    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

Пример 2

program Trap_2
! решение системы уравнений 2-го порядка
! y'' = y' + 2y - 4z•exp(-2x) - 1
! z'' = 2z' + (y - x)•exp(3x)
! с начальными условиями
! y(0)=1, y'(0)=0, z(0)=0, z'(0)=0.5
! методом трапеций с постоянным шагом
use NML
implicit none

    integer, parameter:: m=4, N=2048, NW=1
    character:: met
    real:: Eps
    real:: A, B, H, Y(m), E(m)
    integer:: NE
    integer:: cnt=0, cnt2=0
    real C0, C1, Lam
    integer i
!begin
    met='T'
    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

Вернуться к оглавлению Скачать DE23

DE26

Программа DE26 решает систему из m обыкновенных дифференциальных уравнений первого порядка y_i'(x) = f_i(x, y₁(x), y₂(x), ..., y_m(x)), i=1, 2, ..., m по формуле дифференцирования назад 2-го поядка аппроксимации:
y_n+2' = (4/3)y_n+1 - (1/3)y_n + (2/3)hf(x_n+2, y_n+2)
Для повышения устойчивости итерационный процесс для нахождения y_n+2(x) выполняется по методу Ньютона. Первый шаг программа выполняет неявным методом Эйлера. На следующих шагах предсказание начального приближения осуществляется вычислением экстраполяционного многочлена Эрмита.
На каждом шаге выполняется от 1 до 3 итераций. Итерации прерываются, если каждая из компонент вектоа решения становится меньше заданного числа Eps. Если после выполнения 3-х итераций какая-либо из компонент не достигла заданной точности, специальная переменная nErr увеличивается на 1. На выходе из программы эта переменная содержит общее число шагов, на которых не была достигнута заданная точность Eps.
Для реализации итераций по методу Ньютона на каждой итерации требуется обращение матрицы Якоби, что для сложных систем уравнений может быть достаточно трудоёмко. Если от итерации к итерации матрица Якоби меняется незначительно, заданием специального значения в переменной nWex можно оставить обращение матрицы только перед первой итерацией, а в остальных итерациях пользоваться уже вычисленным значением. Это позволяет значительно уменьшить объём вычислений. В некоторых случаях (для систем уравнений с постоянными коэффициентами) матрица Якоби не зависит от аргумента и функций, т.е. является числовой матрицей. Тогда матрица Якоби обращается только один раз в начале выполнения программы и далее на всём выполнении программы используется уже вычисленное значение.
Алгоритм абсолютно-устойчив и может применяться для решения жёстких уравнений и систем уравнений с прозвольной величиной шага интегрирования.
В процессе вычислений программа 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
! с начальными условиями y(0)=2, z(0)=18

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

Вернуться к оглавлению Скачать DE26

DE53

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

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

D_nk₁ = h•f(y_n)
D_nk₂ = h•f(y_n + β₂₁k₁)
D_nk₃ = h•f(y_n + β₃₁k₁ + β₃₂k₂)
y_n+1 = y_n + p₁k₁ + p₂k₂ + p₃k₃

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

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

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

Норма ║ξ║ расчитывается по формуле
║ξ║ = max_1≤j≤m{|ξ_j|/(|y_n^j| + r)}
где y_n^j - j-ая компонента решения, r – малый положительный параметр. Если |y_n^j| < 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-ой зависимой функции запоминается в элементе DYJ(i,j);
Real X, Xout - начальная и конечная точки интегрирования, после удачного (не аварийного) завершения прграммы X = Xout;
Real H - начальный шаг интегрирования, с которым программе предлагается начать вычисления. Допустимо задать абсолютное значение, даже если шаг должен быть отрицательным. На выходе переменная содержит величину шага в конце промежутка интегрирования;
Real Y(1:M) - массив, который на входе должен содержать начальные значения функций в точке X, на выходе будет содержать решение системы в точке Xout. Число M задает размерность решаемой системы уравнений;
Real Eps - требуемая точность интегрирования на каждом шаге. Так как алгоритм имеет 3-й порядок аппроксимации, оптимальным значением будет ε = 10^-4;
Integer nWex - задаёт режим вычисления Якобиана, может принимать значения 0 или 1; 0 - Якобиан вычисляется только один раз в начале программы; 1 - вычисляется на каждом шаге.
Integer Error - индикатор ошибки.
Error=0, если интегрирование системы успешно закончено и решение вычислено с заданной точностью Eps на каждом шаге;
Error=1, если интегрирование нельзя начать, т.к. значение |Xout-X| меньше минимально-допустимой величины шага h_min, которая вычисляется программой. Параметр H приобретает значение h_min;
Error=2, если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины ε_min, которая вычисляется программой. Параметр Eps приобретает значение ε_min;
Error=65, если программа прекратила вычисления, т.к. решение не может быть вычислено с заданной точностью Eps. Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Решение в точке X должно быть вычислено правильно.
В процессе вычислений программа DE53 вызывает программу AIG2R.

Пример 1

program Rosenbrock1
! решение автономной системы уравнений 1-го порядка
! y' = -20y + z
! z' = 19y - 2z
! с начальными условиями y(0)=2, z(0)=18
! методом Розенброка 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
! с начальными значениями
! y1(0)=1, y2(0)=0, y3(0)=1
! методом Розенброка 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

Вернуться к оглавлению Скачать DE53

DDRIV2

Программа DDRIV2 решает систему из n обыкновенных дифференциальных уравнений первого порядка y_i'(x) = f_i(y₁(x), y₂(x), ..., y_m(x)), i=1, 2, ..., n. Режимы программы позволяют решать как жёсткие, так и не жёсткие системы дифференциальных уравнений. Для не жёстких систем применяются методы Адамса от 1 до 6 порядка, для жёстких методы Гира от 1 до 5 порядка. В обоих случаях шаг интегрирования выбирается автоматически на основе оценки погрешности. Программа является самостартующей.
Система решается итерациями по методу Ньютона. Якобиан может быть задан в соответствующей подпрограмме или вычисляться с помощью конечных разностей. Матрица Якоби может иметь как обычную квадратную, так и ленточную структуру.

В процессе вычислений DDRIV2 использует большое количество вспомогательных подпрограмм, все подпрограммы объеденены в модуле DDRIV, отличным от модуля NML.
Все вещественные переменные в модуле объявлены как переменные двойной точности DOUBLE PRECISION

Вызов программы

CALL DDRIV2 (N, T, Y, F, TOUT, MSTATE, NROOT, EPS, EWT, & MINT, MITER, WORK, LENW, IWORK, LENIW, G,JACOBN,FA,USERS)

Параметры программы

N - (На входе) Число дифференциальных уравнений.
T - Независимая переменная. На входе T задаёт начальную точку интегрирования, на выходе - точка, в которой найдено решение системы.
Y - Массив размерности не меньше N, который на входе должен содержать начальные значения функций в точке T, на выходе будет содержать решение системы в точке TOUT. Массив Y передаётся как параметр в подпрограмму F и в подпрограмму-функцию G. Информация, необходимая для F и G, может быть размещена в этом массиве начиная с (N+1)-го элемента. (Замечание. Изменения, вносимые пользователем в первые N элементов массива, проявятся лишь в случае повторного обращения к программе, т.е. после изменения MSTATE на +1 или -1).
F - Подпрограмма пользователя для вычисления правых частей системы дифференциальных уравнений YDOT = F(Y,T). Здесь Y - массив длины не меньше N, YDOT - массив длины N. Подпрограмма должна иметь вид:

	subroutine F(N, T, Y, YDOT)
	integer N
	double precision T, Y(*), YDOT(*)
	  YDOT(1) = <функция № 1 от T, Y(1), Y(2), ..., Y(N)>
		. . . . . . . . . . . . . . . . . . . .
		. . . . . . . . . . . . . . . . . . . .
	  YDOT(N) = <функция № N от T, Y(1), Y(2), ..., Y(N)>
	end subroutine F

Если пользователь хочет прервать выполнение программы, ему нужно в подпрограмме F присвоить переменной N значение 0, DDRIV2 вернёт управление в программу, вызвавшую DDRIV2, при этом параметр MSTATE примет значение +6(-6). Смена значения N в F не повлияет на значение параметра N, указанного в списке параметров DDRIV2.
TOUT - (На входе) Точка, в которой требуется вычислить решение.
MSTATE - Указатель целого типа, описывающий режим интегрирования. Пользователь может задать MSTATE равным +1 или -1. Если MSTATE положителен, то программа, миновав при интерполировании точку TOUT, проинтерполирует решение. Это наиболее эффективный режим. Если MSTATE отрицателен, то программа будет подбирать шаги интегрирования, чтобы попасть точно в точку TOUT (режим полезен в случае сингулярности правее точки TOUT). Значения MSTATE могут иметь следующий смысл:

(На входе) Такое значение должно быть задано пользователем при первом обращении к программе. При всех последующих обращениях пользователь должен следить за величиной MSTATE, возвращаемой программой. До следующего обращения к DDRIV2 можно только менять знак MSTATE. (Для удобства пользователя, который хотел бы отказаться от задания начальных условий, DDRIV2 допускает вызов со значениями параметров MSTATE = +1(-1) и TOUT = T. В таком случае подпрограмма вернёт MSTATE без изменений, т.е. MSTATE = +1(-1).
(На выходе) Успешное завершение интегрирования. Если требуется продолжить интегрирование (в том же направлении), нужно просто увеличить TOUT и вызвать программу снова.
(На выходе) Неудачное завершение. Выполнено 1000 шагов интегрирования, но точка TOUT не достигнута. Пользователь может продолжить интегрирование, снова обратившись к DDRIV2. Если это не ошибка в постановке задачи, то наиболее вероятной причиной может быть попытка решить жёсткую систему уравнений в режиме нежёсткого интегрирования. (См. ниже описание параметра MINT.)
(На выходе) Неудачное завершение. Затребована слишком высокая точность. EPS было увеличено до приемлемого уровня, который рассчитала программа. Пользователь может продолжить интегрирование, снова обратившись к DDRIV2.
(На выходе) Корень был найден в точке, лежащей левее точки TOUT. Пользователь может продолжить интегрирование к точке TOUT, снова обратившись к DDRIV2.
(На выходе) Неудачное завершение. В подпрограмме F для N было установлено значение 0.
(На выходе) Неудачное завершение. В подпрограмме-функции 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 и Якоби являются ленточными.
WORK
LENW - (На входе) WORK вещественный массив длины LENW, используемый внутри программы в качестве рабочего. Пользователь должен отвести в вызывающей программе под этот массив память посредством описания
REAL WORK(...) или REAL(kind=8) WORK(...)
Длина массива WORK должна быть не меньше, чем
16 × T + 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(...)
Длина IWORK должна быть не меньше, чем
21, если MINT = 1,
N + 21, если MINT = 2 или 3
и установленный размер памяти следует присвоить параметру LENIW. Содержимое массива IWORK не должно изменяться между обращениями к DDRIV2.
G - Вещественная подпрограмма-функция, которую пользователь должен написать, если NROOT не есть 0. Функция 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 (2)
	      G= ...
	    .................
	    .................
	    case default
	      G= ...
	  end select
	end function G

Здесь Y - массив длины не меньшей N; его первые N компонент содержат решение в точке T. Пользователь не должен изменять эти значения. Фактический размер массива Y определяется по его описанию в программе пользователя, вызывающей DDRIV2. Однако, если пользователь желает прервать вычисления, т.е. передать управление в программу, вызвавшую DDRIV2, ему следует присвоить N нулевое значение. DDRIV2 сигнализирует об этом, возвратив MSTATE со значением, равным 7(-7). В этом случае номер обработанного в G уравнения будет размещён в шестой компоненте массива IWORK. Смена значения N в G не повлияет на значение параметра N, указанного в списке параметров DDRIV2.
JACOBN - Подпрограмма, предоставленная пользователем, если значение MITER равно 1 или 4. Подпрограмма вычисляет частные производные I-й правой части по отношению к J-й зависимой переменной. Для MITER = 1 матрица Якоби квадратная размером N×N, для MITER = 4 матрица ленточная с ML лентами ниже и MU лентами выше главной диагонали. Более подробно смотрите описание формальных параметров к программе DDRIV3.
FA - Подпрограмма, предоставленная пользователем, если значение IMPL равно 1 или 2, а значение MITER не равно 3. Эта подпрограмма вычисляет массив A, где A*dY(I)/dT = F(Y(I),T). Более подробно смотрите описание формальных параметров к программе DDRIV3.
USERS - Подпрограмма, предоставляемая пользователем, если MITER = 3. Подпрограмма USERS вызывается DDRIV3, когда необходимо решить определенные линейные задачи. Более подробно смотрите описание формальных параметров к программе DDRIV3.

Дополнительная информация для пользователя

Первые три элемента массива WORK и первые пять элементов массива IWORK будет помещена следующая статистическая информация:
AVGH - среднее значение использованных шагов,
HUSED - последний (успешный) использованный шаг,
AVGORD - среднее значение использованных порядков,
IMXERR - индекс того элемента вектора решения, на котором достигается максимальная ошибка при последней проверке,
NQUSED - последний (успешно) использованный порядок,
NSTEP - число шагов, выполненных после последней инициализации,
NFE - число вычислений правых частей,
NJE - число вычислений матриц-якобианов.
В процессе вычислений DDRIV2 вызывает подпрограммы DDRIV3 и XERROR.

В связи с тем, что пакет содержит большое количество программ и имеет большой объём, пакет размещён в отельном модуле DDRIV и не использует программ из модуля NML. Тексты программ переведены с Фотран 77 на Фортран 90. Изменения незначительные. В число формальных параметров добавлена переменная MITER, теперь пользователь должен сам указать при вызове правильное значение этой переменной. Подпрограммы d1mach и i1mach были изменены для использования встроенных функций, доступных в Fortran 90. Данные программы представлены на Международном семинаре «Современные направления в численном программном обеспечении и в высокопроизводительных вычислениях», 19–20 октября 1995 г., Киото, Япония.
Исходные тексты программ пакета на языке Фортран 77 предоставлены автору Василием Рябовым <****@****>, за что автор выражает ему свою искреннюю благодарность. Краткое описание программы DDRIV2 из пакета можно найти в [Б19].

Пример 1

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             ! Вычисление точных значений в точке 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

Пример 2

program TestDDRIV2_G
    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             ! Вычисление точных значений в точке 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

Вернуться к оглавлению Скачать DDRIV

DT10

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

Вызов программы

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

Параметры программы

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

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

Пример

program test_dt10
! Краевая задача для уравнения
! y" - 2x·y' - 2·y = 4x
! y(0) = 1;  y(1) = ℮ - 1
use NML
implicit none

integer, parameter:: n=21
real X(n), Y(n), E(n)
real A, B, YA, YB
real h, R
integer i
!begin
    A=0.0; B=1.0
    ! начальные значения
    YA=exp(A*A)-A; YB=exp(B*B)-B

    ! точные значения
    h=(B-A)/real(n-1)
    do i=1, n
        R=A+h*real(i-1); E(i)=exp(R*R)-R
    end do

    ! X(1)=A;  X(n)=B
    ! решение задачи методом прогонки
    call DT10R(Q, P, F, A, B, YA, YB, X, Y)

    print 10
    print 11, (X(i), Y(i), E(i), i=1,n,2)
10  format(/10X,'X',14X,'Y',17X,'E'/)
11  format(4X,F10.6,4X,E14.7,4X,E14.7)
stop

contains

real function Q(X)
    real, intent(in):: X
    Q=-2.0*X
end function Q

real function P(X)
    real, intent(in):: X
    P=-2.0
end function P

real function F(X)
    real, intent(in):: X
    F=4.0*X
end function F

end program test_dt10

          X              Y                 E

      0.000000     0.1000000E+01     0.1000000E+01
      0.100000     0.9101487E+00     0.9100502E+00
      0.200000     0.8409864E+00     0.8408108E+00
      0.300000     0.7944074E+00     0.7941743E+00
      0.400000     0.7737821E+00     0.7735109E+00
      0.500000     0.7843144E+00     0.7840254E+00
      0.600000     0.8336133E+00     0.8333294E+00
      0.700000     0.9325696E+00     0.9323162E+00
      0.800000     0.1096676E+01     0.1096481E+01
      0.900000     0.1348016E+01     0.1347908E+01
      1.000000     0.1718282E+01     0.1718282E+01

Вернуться к оглавлению Скачать DT10

DT11

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

Вызов программы

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

Параметры программы

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

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

Пример 1

program test_dt11
! Краевая задача для уравнения
! y" + 2x/(x² + 1)·y' + 2/(x² + 1)·y = 4x/(x² + 1)³
! y'(-1) - y(-1) = 0
! y'(1) + y(1) = 0
use NML
implicit none

integer, parameter:: n=81
real X(n), Y(n), E(n)
real A, B, Alpha, Alpha1, Beta, Beta1
real h, R
integer i
!begin
    A=-1.0; B=1.0
    ! начальные значения
    Alpha=-1.0; Alpha1=0.0
    Beta=1.0; Beta1=0.0

    ! точные значения
    h=(B-A)/real(n-1)
    do i=1, n
        R=A+h*real(i-1); E(i)=1.0/(R*R+1)
    end do

    ! X(1)=A;  X(n)=B
    ! решение задачи методом прогонки
    call DT11R(Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1, X, Y)

    print 10
    print 11, (X(i), Y(i), E(i), i=1,n,8)
10  format(/10X,'X',14X,'Y',17X,'E'/)
11  format(4X,F10.6,4X,E14.7,4X,E14.7)
stop

contains

real function Q(X)
    real, intent(in):: X
    Q=2.0*X/(X*X+1)
end function Q

real function P(X)
    real, intent(in):: X
    P=2.0/(X*X+1)
end function P

real function F(X)
    real, intent(in):: X
    F=4.0*X*X/(X*X+1)**3
end function F

end program test_dt11

          X              Y                 E

     -1.000000     0.4998799E+00     0.5000000E+00
     -0.800000     0.6096135E+00     0.6097562E+00
     -0.600000     0.7351372E+00     0.7352942E+00
     -0.400000     0.8619106E+00     0.8620690E+00
     -0.200000     0.9613886E+00     0.9615385E+00
      0.000000     0.9998574E+00     0.1000000E+01
      0.200000     0.9613931E+00     0.9615385E+00
      0.400000     0.8619165E+00     0.8620689E+00
      0.600000     0.7351429E+00     0.7352940E+00
      0.800000     0.6096182E+00     0.6097561E+00
      1.000000     0.4998826E+00     0.5000000E+00

Пример 2

program boundary_task_2
! Краевая задача для уравнения
! y" - 2x·y' - 2·y = 4x
! y'(0) + 3·y(0) = 2
! y'(1) + (4℮ - 3)/(℮ - 1)·y(1) = 6℮ - 4
use NML
implicit none
integer, parameter:: n=20
real X(0:n), Y(0:n), E(0:n)
! размерность массивов X и Y 21 элемент
real A, B, Alpha, Alpha1, Beta, Beta1
real h, R, e1
integer i
!begin
    A=0.0; B=1.0
    ! начальные значения
    Alpha=3.0; Alpha1=2.0
    e1=exp(1.0)
    Beta=(4.0*e1-3.0)/(e1-1.0); Beta1=6.0*e1-4.0

    ! точные значения
    h=(B-A)/real(n)
    do i=0, n
        R=A+h*real(i); E(i)=exp(R*R)-R
    end do

    ! X(0)=A;  X(n)=B
    ! решение задачи методом прогонки
    call DT11R(Q, P, F, A, B, Alpha, Alpha1, Beta, Beta1, X, Y)

    print 10
    print 11, (X(i), Y(i), E(i), i=0,n,2)
10  format(/10X,'X',14X,'Y',17X,'E'/)
11  format(4X,F10.6,4X,E14.7,4X,E14.7)
stop

contains

real function Q(X)
    real, intent(in):: X
    Q=-2.0*X
end function Q

real function P(X)
    real, intent(in):: X
    P=-2.0
end function P

real function F(X)
    real, intent(in):: X
    F=4.0*X
end function F

end program boundary_task_2

          X              Y                 E

      0.000000     0.1000236E+01     0.1000000E+01
      0.100000     0.9102049E+00     0.9100502E+00
      0.200000     0.8408605E+00     0.8408108E+00
      0.300000     0.7940893E+00     0.7941743E+00
      0.400000     0.7732548E+00     0.7735109E+00
      0.500000     0.7835505E+00     0.7840254E+00
      0.600000     0.8325718E+00     0.8333294E+00
      0.700000     0.9311924E+00     0.9323162E+00
      0.800000     0.1094882E+01     0.1096481E+01
      0.900000     0.1345695E+01     0.1347908E+01
      1.000000     0.1715279E+01     0.1718282E+01

Вернуться к оглавлению Скачать DT11