Автор cайта: Владимир Потемкин fortran-90@yandex.ru |
|
Численное интегрирование
DF40. Вычисление множества значений интегралов для функции, заданной таблично в неравноотстоящих точках.
DF42. Вычисление множества значений интегралов для функции, заданной таблично в равноотстоящих точках по трем точкам.
DF43. Вычисление множества значений интегралов для функции, заданной таблично в равноотстоящих точках по пяти точкам.
DF47. Вычисление множества значений интегралов для функции, заданной значениями самой функции и значениями ее первой прозводной в неравноотстоящих точках.
DF48. Вычисление множества значений интегралов для функции, заданной значениями самой функции и значениями ее первой прозводной в равноотстоящих точках.
DF60. Вычисление значения интеграла для функции, заданной аналитически по формуле Боде с автоматическим выбором величины шага интегрирования.
DF63. Вычисление значения интеграла для функции, заданной таблично в равноотстоящих точках с шагом h по формуле Симпсона.
DF64. Вычисление значения интеграла для функции, заданной таблично в равноотстоящих точках с шагом h по формуле Грегори пятого порядка точности.
GA2 - GA64. Вычисление значения интеграла для функции, заданной аналитически по формулам Гаусса с числом узлов от 2 до 64.
DF67. Вычисление значения интеграла для функции, заданной таблично в неравноотстоящих точках по формуле трапеций.
DF40
Программа DF40 вычисляет множество z1, z2, ..., zn значений интегралов
xi
zi = ∫ y(x)dx (i=1, 2, ..., n)
x1
от функции y(x), заданной множествами неравноотстоящих точек x1, x2, ..., xn упорядоченных по возрастанию или убыванию значений аргумента и y1, y2, ..., yn соответствующих значений функции.
Учитывая, что z1=0, последовательные значения интегралов вычисляются по формуле трапеций:
zi = zi-1 + ½(xi-xi-1)(yi+yi-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 вычисляет множество z1, z2, ..., zn значений интегралов
xi
zi = ∫ y(x)dx (i=1, 2, ..., n)
x1
от функции y(x), заданной множеством y1, y2, ..., yn ее значений в равноотстоящих точках x1, x2, ..., xn с шагом изменения аргумента h. Предполагается, что функция непрерывна и, по крайней мере, дважды дифференцируема.
Учитывая, что z1=0, последовательные значения интегралов вычисляются по формуле трапеций:
zi = zi-1 + ½ h(yi+yi-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 вычисляет множество z1, z2, ..., zn значений интегралов
xi
zi = ∫ y(x)dx (i=1, 2, ..., n)
x1
от функции y(x), заданной множеством y1, y2, ..., yn ее значений в равноотстоящих точках x1, x2, ..., xn с шагом изменения аргумента h. Предполагается, что функция непрерывна и дифференцируема пять раз.
В зависимости от значения n для вычисления интегралов при различных значениях i в программе используется формула трапеций (1), формула Симпсона (2) и (3), формула "трех восьмых" (4), (5) и (6) или формула Боде (7), (8), (9) и (10).
z1 = 0 | |
z2 = (h/2)•(y1+y2) | (1) |
z2 = (h/12)•(5y1+8y2-y3) | (2) |
z3 = (h/3)•(y1+4y2+y3) | (3) |
z2 = (h/24)•(9y1+19y2-5y3+y4) | (4) |
z3 = (h/3)•(y1+4y2+y3) | (5) |
z4 = (3h/8)•(y1+3y2+3y3+y4) | (6) |
z2 = (h/720)•(251y1+646y2-264y3+106y4-19y5) | (7) |
z3 = (h/90)•(29y1+124y2+24y3+4y4-y5) | (8) |
z4 = (h/80)•(27y1+102y2+72y3+42y4-3y5) | (9) |
zi = zi-4 + (2h/45)•(7yi-4+32yi-3+12yi-2+32yi-1+7yi) (i=5, ..., n) | (10) |
Вызов программы
CALL DF43(Y, H, IY)
Параметры программы
Y, H - входные параметры;
IY - выходной параметр;
Real Y(1:N) - массив, содержащий значения функции;
Real H - шаг изменения аргумента;
Real IY(1:N) - массив, содержащий вычисленные значения интеграла. При вызове массив можно совмещать в памяти с массивом Y.
Вернуться к оглавлению Скачать DF43
DF47
Программа DF47 вычисляет множество z1, z2, ..., zn значений интегралов
xi
zi = ∫ y(x)dx (i=1, 2, ..., n)
x1
от функции y(x), заданной множеством x1, x2, ..., xn упорядоченных по возрастанию или убыванию значений аргумента и множествами y1, y2, ..., yn и y'1, y'2, ..., y'n соответствующих значений функции и ее первой производной. Предполагается, что интегрируемая функция непрерывна и дифференцируема четыре раза.
Учитывая, что z1=0, последовательные значения интегралов вычисляются по формуле:
zi = zi-1 + (1/2)(xi-xi-1)[yi-1 + yi + (1/6)(xi-xi-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 вычисляет множество z1, z2, ..., zn значений интегралов
xi
zi = ∫ y(x)dx (i=1, 2, ..., n)
x1
от функции y(x), заданной множеством y1, y2, ..., yn ее значений и множеством y'1, y'2, ..., y'n значений ее производной в равноотстоящих точках с шагом изменения аргумента h. Предполагается, что интегрируемая функция непрерывна и дифференцируема четыре раза.
Учитывая, что z1=0, последовательные значения интегралов вычисляются по формуле:
zi = zi-1 + (1/2)h[yi-1 + yi + (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), заданную аналитически с автоматическим выбором величины шага исходя из оценки погрешности интегрирования на каждом шаге. Погрешность интегрирования оценивается на основе анализа двух квадратурных формул
I0 = (h/90)[7f(x) + 32f(x+h/4) + 12f(x+h/2) + 32f(x+3h/4) + 7f(x+h)]
I1 = (h/9450)[455f(x) + 1536f(x+h/8) + 3584f(x+3h/8) - 420f(x+h/2) +
+ 3584f(x+3h/4) + 711f(x+h)],
для которых справедливы следующие асимптотические оценки погрешности:
R0 = I - I0 = -(8/945)•(h/4)7•f(VI)(ξ)
R1 = I - I1 = -(26/4725)•(h/4)7•f(VI)(ξ),
где I обозначает точное (неизвестное) значение интеграла. За приближение к точному значению интеграла принимается I0, а общая погрешность интегрирования оценивается по правилу Рунге
I1 - I0 = R0 - R1 = [1 - (26/4725)•(945/8)]•[-(8/945)•(h/4)7•f(VI)(ξ)] = (7/20)R0,
откуда R0 = (20/7)(I1 - I0), где
I1 - I0 = (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)].
Шаг интегрирования программа выбирает автоматически. Если вычисленная погрешность R0 превосходит заданное значение ε, то шаг интегрирования уменьшается вдвое, причем вычисленные значения функции запоминаются для последующего использования. Если погрешность R0 удовлетворяет заданной точности, то вычисленное значение I0 принимается за значение интеграла на данном шаге, а величина следующего шага hnew рассчитывается по формуле:
7 ______
hnew = hold•√ ε/R0.
Программа наиболее эффективна для быстро меняющихся (быстро растущих или убывающих) или осциллирующих функций.
Вызов программы
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| меньше минимально-допустимой величины шага hmin, которая вычисляется программой. Параметр H приобретает значение hmin;
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)•∑ ½ck•f(½(b-a)dk + ½(b+a)),
k=1
где dk - абсциссы (узлы) квадратурной формулы Гаусса n-го порядка для промежутка [-1,1], ck - коэффициенты этой формулы.
Число узлов квадратурной формулы, используемых для вычисления интеграла, отражено в названии программы.
Вызов программы
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] на неравномерной сетке xi, i = 1, ..., N, по формуле трапеций.
A = x1 и B = xN - нижний и верхний пределы интегрирования.
Вызов программы
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
|