Численное интегрирование. Библиотека численных методов на языке Fortran 90

DF40

Программа DF40 вычисляет множество z₁, z₂, ..., z_n значений интегралов
_{x_i}
z_i = ∫ y(x)dx (i=1, 2, ..., n)
^x₁
от функции y(x), заданной множествами неравноотстоящих точек x₁, x₂, ..., x_n упорядоченных по возрастанию или убыванию значений аргумента и y₁, y₂, ..., y_n соответствующих значений функции.
Учитывая, что z₁=0, последовательные значения интегралов вычисляются по формуле трапеций:
z_i = z_i-1 + ½(x_i-x_i-1)(y_i+y_i-1) (i=2, 3, ..., n).

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

CALL DF40(X, Y, IY)

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

X, Y - входные параметры;
IY - выходной параметр;

Real X(1:N), Y(1:N) - массивы, задающие значения аргумента и функции.
Real IY(1:N) - массив, содержащий вычисленные значения интеграла. При вызове массив можно совмещать в памяти с массивом X или Y.

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

DF42

Программа DF42 вычисляет множество z₁, z₂, ..., z_n значений интегралов
_{x_i}
z_i = ∫ y(x)dx (i=1, 2, ..., n)
^x₁
от функции y(x), заданной множеством y₁, y₂, ..., y_n ее значений в равноотстоящих точках x₁, x₂, ..., x_n с шагом изменения аргумента h. Предполагается, что функция непрерывна и, по крайней мере, дважды дифференцируема.
Учитывая, что z₁=0, последовательные значения интегралов вычисляются по формуле трапеций:
z_i = z_i-1 + ½ h(y_i+y_i-1) (i=2, 3, ..., n).

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

CALL DF42(Y, H, IY)

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

Y, H - входные параметры;
IY - выходной параметр;

Real Y(1:N) - массив, содержащий значения функции;
Real H - шаг изменения аргумента;
Real IY(1:N) - массив, содержащий вычисленные значения интеграла. При вызове массив можно совмещать в памяти с массивом Y.

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

DF43

Программа DF43 вычисляет множество z₁, z₂, ..., z_n значений интегралов
_{x_i}
z_i = ∫ y(x)dx (i=1, 2, ..., n)
^x₁
от функции y(x), заданной множеством y₁, y₂, ..., y_n ее значений в равноотстоящих точках x₁, x₂, ..., x_n с шагом изменения аргумента h. Предполагается, что функция непрерывна и дифференцируема пять раз.
В зависимости от значения n для вычисления интегралов при различных значениях i в программе используется формула трапеций (1), формула Симпсона (2) и (3), формула "трех восьмых" (4), (5) и (6) или формула Боде (7), (8), (9) и (10).

`z₁ = 0`
`z₂ = (h/2)•(y₁+y₂)`	(1)
`z₂ = (h/12)•(5y₁+8y₂-y₃)`	(2)
`z₃ = (h/3)•(y₁+4y₂+y₃)`	(3)
`z₂ = (h/24)•(9y₁+19y₂-5y₃+y₄)`	(4)
`z₃ = (h/3)•(y₁+4y₂+y₃)`	(5)
`z₄ = (3h/8)•(y₁+3y₂+3y₃+y₄)`	(6)
`z₂ = (h/720)•(251y₁+646y₂-264y₃+106y₄-19y₅)`	(7)
`z₃ = (h/90)•(29y₁+124y₂+24y₃+4y₄-y₅)`	(8)
`z₄ = (h/80)•(27y₁+102y₂+72y₃+42y₄-3y₅)`	(9)
`z_i = z_i-4 + (2h/45)•(7y_i-4+32y_i-3+12y_i-2+32y_i-1+7y_i) (i=5, ..., n)`	(10)

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

CALL DF43(Y, H, IY)

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

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

DF47

Программа DF47 вычисляет множество z₁, z₂, ..., z_n значений интегралов
_{x_i}
z_i = ∫ y(x)dx (i=1, 2, ..., n)
^x₁
от функции y(x), заданной множеством x₁, x₂, ..., x_n упорядоченных по возрастанию или убыванию значений аргумента и множествами y₁, y₂, ..., y_n и y'₁, y'₂, ..., y'_n соответствующих значений функции и ее первой производной. Предполагается, что интегрируемая функция непрерывна и дифференцируема четыре раза.
Учитывая, что z₁=0, последовательные значения интегралов вычисляются по формуле:
z_i = z_i-1 + (1/2)(x_i-x_i-1)[y_i-1 + y_i + (1/6)(x_i-x_i-1)(y'_i-1-y'_i)] (i=2, 3, ..., n).

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

CALL DF47(X, Y, DY, IY)

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

X, Y, DY - входные параметры;
IY - выходной параметр;

Real X(1:N), Y(1:N), DY(1:N) - массивы, содержащие значения аргумента, функции и ее первой производной;
Real IY(1:N) - массив, содержащий вычисленные значения интеграла. При вызове массив можно совмещать в памяти с массивом X, Y или DY.

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

DF48

Программа DF48 вычисляет множество z₁, z₂, ..., z_n значений интегралов
_{x_i}
z_i = ∫ y(x)dx (i=1, 2, ..., n)
^x₁
от функции y(x), заданной множеством y₁, y₂, ..., y_n ее значений и множеством y'₁, y'₂, ..., y'_n значений ее производной в равноотстоящих точках с шагом изменения аргумента h. Предполагается, что интегрируемая функция непрерывна и дифференцируема четыре раза.
Учитывая, что z₁=0, последовательные значения интегралов вычисляются по формуле:
z_i = z_i-1 + (1/2)h[y_i-1 + y_i + (1/6)h(y'_i-1-y'_i)] (i=2, 3, ..., n).

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

CALL DF48(Y, DY, H, IY)

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

Y, DY, H - входные параметры;
IY - выходной параметр;
Real Y(1:N), DY(1:N) - массивы, содержащие значения функции и ее первой производной;;
Real H - шаг изменения аргумента;
Real IY(1:N) - массив, содержащий вычисленные значения интеграла. При вызове массив можно совмещать в памяти с массивом Y или DY.

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

DF60

Программа DF60 интегрирует функцию f(x), заданную аналитически с автоматическим выбором величины шага исходя из оценки погрешности интегрирования на каждом шаге. Погрешность интегрирования оценивается на основе анализа двух квадратурных формул

I₀ = (h/90)[7f(x) + 32f(x+h/4) + 12f(x+h/2) + 32f(x+3h/4) + 7f(x+h)]
I₁ = (h/9450)[455f(x) + 1536f(x+h/8) + 3584f(x+3h/8) - 420f(x+h/2) +
+ 3584f(x+3h/4) + 711f(x+h)],

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

R₀ = I - I₀ = -(8/945)•(h/4)⁷•f^(VI)(ξ)
R₁ = I - I₁ = -(26/4725)•(h/4)⁷•f^(VI)(ξ),

где I обозначает точное (неизвестное) значение интеграла. За приближение к точному значению интеграла принимается I₀, а общая погрешность интегрирования оценивается по правилу Рунге

I₁ - I₀ = R₀ - R₁ = [1 - (26/4725)•(945/8)]•[-(8/945)•(h/4)⁷•f^(VI)(ξ)] = (7/20)R₀,

откуда R₀ = (20/7)(I₁ - I₀), где

I₁ - I₀ = (h/4725)[-140f(x) + 768f(x+h/8) + 1792f(x+3h/8) - 1680f(x+h/4) -
- 840f(x+h/2) + 112f(x+3h/4) - 12f(x+h)].

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

Программа наиболее эффективна для быстро меняющихся (быстро растущих или убывающих) или осциллирующих функций.

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

CALL DF60(Fun, X, Xout, H, Eps, SI, SE, Error)

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

Fun, Xout - входные параметры;
X, H, Eps - входные и выходные параметры;
SI, SE, Error - выходные параметры;

Real Fun(Real X) - интегрируемая функция f(x);
Real X, Xout - начальная и конечная точки интегрирования;
Real H - начальный шаг интегрирования, с которым программе предлагается начать вычисления. Допустимо задать абсолютное значение, даже если шаг должен быть отрицательным. На выходе переменная содержит величину шага в конце промежутка интегрирования;
Real Eps - требуемая точность интегрирования на каждом шаге;
Real SI - вычисленное значение интеграла;
Real SE - сумма оценок погрешностей интегрирования на каждом шаге;
Integer Error - индикатор ошибки;
Error=0, если интеграл вычислен и погрешность интегрирования на каждом шаге удовлетворяет заданной точности;
Error=1, если интегрирование нельзя начать, т.к. значение |Xout-X| меньше минимально-допустимой величины шага h_min, которая вычисляется программой. Параметр H приобретает значение h_min;
Error=2, если интегрирование нельзя начать, т.к. заданное значение Eps меньше минимально-допустимой величины ε_min, которая вычисляется программой. Параметр Eps приобретает значение ε_min;
Error=65, если программа прекратила вычисления, т.к. интеграл не может быть вычислен с заданной точностью Eps. Переменная X получает значение аргумента, при котором произошло прекращение вычислений. Интеграл от начального значения X до его значения при выходе из программы вычислен правильно.

Пример

! Интерирование функции Fun(x) с автоматическим
! выбором величины шага интегрирования.

program Test_DF60
use NML
implicit none
real:: A, B, H, E, Eps, S1, S2
integer Er, Count
!begin
  Count=0
  Eps=1.0E-7
  A=-1.0; B=1.0; H=0.0625
  E=10.0*(atan(10.0)-atan(-10.0))
  call DF60(Fun, A, B, H, Eps, S1, S2, Er)
  print 10, Er, E, S1, S2, Count
10 format(/'      Error =',I3/'    E =',F12.7    &
          /'   S1 =',F12.7/'   S2 =',E16.7/'      Count =',I5)

contains

real function Fun(x)
real, intent(in):: x
  Fun=1.0/(x*x+0.01)
  Count=Count+1
  return
end function Fun

end program Test_DF60

      Error =  0
    E =  29.4225521
   S1 =  29.4225521
   S2 =   0.8375608E-05
      Count =  121

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

DF63

Программа DF63 интегрирует функцию f(x), заданную таблично массивом значений F[1..N] на равномерной сетке x = A + (i - 1)*H, i = 1, ..., N с шагом H по квадратурной формуле Симпсона, точной для многочленов второй степени.
A - нижний предел интегрирования, верхний предел интегрирования определяется из равенства B = A + H*(N - 1)

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

DF63(F, H)

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

F, H - входные параметры;
Real DF63 - возвращаемое значение, содержащее вычисленное значение интеграла;
Real F(1:N) - значения интегрируемой функции;
Real H - шаг сетки.

Пример

! Интегрирование функции F, заданнной таблично на равномерной сетке, по формуле Симпсона
! A, B - нижний и верхний пределы интегрирования

program Test_DF63
    use NML
    implicit none
    integer, parameter:: N=120
    real F(N), A, B, H, X, RF, EF
    data A/-1.0/, B/1.0/
    integer i
!begin
    H=(B-A)/(N-1)
    do i=1, N
        X=A+H*(i-1)
        F(i)=1.0/(X*X+0.01)
    enddo
    RF=DF63(F, H)
    EF=10.0*(atan(10.0)-atan(-10.0))
    print*, 'N=',N
    print*, 'Calculated value',RF
    print*, 'Exact value',EF
!end
end program Test_DF63

Вывод программы:
 N=         40
 Calculated value   29.42225
 Exact value   29.42255

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

DF64

Программа DF64 интегрирует функцию f(x), заданную таблично массивом значений F[1..N] на равномерной сетке x = A + (i - 1)*H, i = 1, ..., N с шагом H, по формуле Грегори пятого порядка точности.
A - нижний предел интегрирования, верхний предел интегрирования определяется из равенства B = A + H*(N - 1)

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

DF64(F, H)

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

F, H - входные параметры;
Real DF64 - возвращаемое значение, содержащее вычисленное значение интеграла;
Real F(1:N) - значения интегрируемой функции;
Real H - шаг сетки.

Пример

! Интегрирование функции F, заданнной таблично, по формуле Грегори
! A, B - нижний и верхний пределы интегрирования

program Test_DF64
    use NML
    implicit none
    integer, parameter:: N=40
    real F(N), A, B, H, X, RF, EF
    data A/-1.0/, B/1.0/
    integer i
!begin
    H=(B-A)/(N-1)
    do i=1, N
        X=A+H*(i-1)
        F(i)=1.0/(X*X+0.01)
    enddo
    RF=DF64(F, H)
    EF=10.0*(atan(10.0)-atan(-10.0))
    print*, 'N=',N
    print*, 'Calculated value',RF
    print*, 'Exact value',EF
!end
end program Test_DF64

Вывод программы:
 N=         40
 Calculated value   29.42224
 Exact value   29.42255

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

GA2 GA4 GA8 GA16 GA32 GA64

Программы интегрируют функцию f(x), заданную аналитически, по квадратурным формулам Гаусса с n узлами в промежутке [a, b]:
_n
I = (b-a)•∑ ½c_k•f(½(b-a)d_k + ½(b+a)),
^k=1

где d_k - абсциссы (узлы) квадратурной формулы Гаусса n-го порядка для промежутка [-1,1], c_k - коэффициенты этой формулы.
Число узлов квадратурной формулы, используемых для вычисления интеграла, отражено в названии программы.

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

GAn(Fun, AX, BX), n равно 2, 4, 8, 16, 32 или 64;

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

Fun, AX, BX - входные параметры;
Real GAn - возвращаемое значение, содержащее вычисленное значение интеграла;
Real Fun(Real X) - интегрируемая функция f(x);
Real AX, BX - начальная и конечная точки интегрирования.

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

DF67

Программа DF67 интегрирует функцию f(x), заданную таблично массивом значений F[1..N] на неравномерной сетке x_i, i = 1, ..., N, по формуле трапеций.
A = x₁ и B = x_N - нижний и верхний пределы интегрирования.

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

DF67(X, F)

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

X, F - входные параметры;
Real DF67 - возвращаемое значение, содержащее вычисленное значение интеграла;
Real X(1:N) - значения аргумента в узлах сетки;
Real F(1:N) - значения подъинтегральной функции;

Пример

! Интегрирование функции F, заданнной таблично на неравномерной сетке, по формуле трапеций
! A, B - нижний и верхний пределы интегрирования

program Test_DF67
    use NML
    implicit none
    integer, parameter:: N=40
    real X(N), F(N), A, B, H, XX, RF, EF, D
    data A/-1.0/, B/1.0/, D/0.125/
    integer i
!begin
    H=(B-A)/(N-1);  D=D*H
    do i=2, N-1
        X(i)=A+H*(i-1)
        if(mod(i,2)==0) then; X(i)=X(i)+D  ! смещение чётных точек вправо
        else; X(i)=X(i)-D  ! смещение нечётных точек влево
        endif
        F(i)=1.0/(X(i)*X(i)+0.01)
    enddo
    X(1)=A;  F(1)=1.0/(X(1)*X(1)+0.01)
    X(N)=B;  F(N)=1.0/(X(N)*X(N)+0.01)
    RF=DF67(X, F)
    EF=10.0*(atan(10.0)-atan(-10.0))
    print*, 'N=',N
    print*, 'Calculated value',RF
    print*, 'Exact value',EF
!end
end program Test_DF67

Вывод программы:
 N=          40
 Calculated value   29.47384
 Exact value   29.42255

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