Zeroin
Программа Zeroin вычисляет корень функции f(x) в заданном интервале [a,b]. Интервал должен содержать не более одного корня функции первой кратности, т.е. f(a) и f(b) должны быть разных знаков.
Для вычисления корня программа использует методы деления отрезка пополам (бисекция), при возмоожности переходит к линейной (метод секущих) или обратной квадратичной интерполяции.
Подробное описание программы "ZEROIN" можно найти в книге [Б3].
Вызов программы
Zeroin(Fun, AX, BX, Eps)
Параметры программы
Fun, AX, BX, Eps - входные параметры;
Real Zeroin - возвращаемое значение;
Real Fun(Real X) - исследуемая функция;
Real AX, BX - границы интервала;
Real Eps - погрешность вычисления корня.
Пример
!Поиск корня нелинейного уравнения
program TestZeroin
use NML
implicit none
integer:: Count=0
real:: A, B, R, Eps
data A/-5.0/, B/5.0/, Eps/1.0E-8/
!begin
R=Zeroin(Fun, A, B, Eps)
print 10, R, Count
10 format(/'Root =',E16.7,' Count =',I3)
contains
real function Fun(x)
real, intent(in):: x
Fun=((x+1.0)*x+2.0)*((x+2.0)*x+2.0)* &
((x-1.0)*x+2.0)*(x-4.0)*((x-8.0)*x+20.0)
Count=Count+1
return
end function Fun
end program TestZeroin
Root = 0.4000000E+01 Count = 13
ZR10
Программа ZR10 вычисляет n корней произвольной функции f(x) при помощи изучения поведения знака функции в интервале [a,b] и дальнейшее их уточнение с точностью до Eps.
Программа "шагает" от точки a к точке b с шагом h. Если функция меняет знак, программа уточняет значение корня с точностью до Eps делением отрезка пополам (бисекция). Используется относительная погрешность (|x - xroot| / |x|) < Eps, в окрестности нуля - абсолютная погрешность |x - xroot| < Eps, а также контролируется значение функции |Fun(x)| < Eps. Не стоит брать значение Eps слишком маленьким. Полученные значения корней следует рассматривать только как приближения к точным значениям и использовать это приближение в других, более быстродействующих, программах (например Zeroin).
b может быть как больше a, так и меньше, в последнем случае программа "шагает" от больших значений аргумента к меньшим. Рекомендуется делать шаг h не больше предполагаемого минимального расстояния между корнями.
Вызов программы
call ZR10(Fun, A, B, H, Eps, C, Cnt)
Параметры программы
Fun, A, B, H, Eps - входные параметры;
Real Fun(Real X) - исследуемая функция;
Real A, B - границы интервала;
Real H - шаг поиска корней;
Real Eps - погрешность вычисления корня;
Real C(N) - массив, в который помещаются найденные корни;
Integer N - количество предполагаемых корней;
Integer Cnt - количество найденных корней.
ZR30
Программа ZR30 вычисляет вещественные корни функции f(x) методом парабол (Мюллера).
На каждой k-ой итерации по трём точкам pk, pk+1, pk+2 строится
интерполяционный полином Лагранжа L(p):
и за следуещее приближение принимается ближайший к pk+2 корень уравнения L(p)=0. Итерационный процесс продолжается, пока приближение к корню не достигнет заданной точности. Если задано начальное приближение к i-му корню ci, то за первые три точки принимаются значения p0=0.9×ci, p1=1.1×ci, p2=1.0×ci. Если ci=0, то в качестве перых трёх точек принимается -1, 1 и 0.
Для облегчения вычислений в программе вводятся следующие переменные:
За новую независимую переменную принимается dk. Тогда корни уравнения L(dk)=0 будут равны
Знак перед фигурной скобкой выбирается так, чтобы модуль знаменателя был наибольшим. Следующая точка в последовательности будет равна
pk+1 = pk + dk+1hk
На каждом шаге вычисляется соотношение | F(pk+1) / F(pk) |. Если это соотношение уменьшается больше чем в 10 раз, то dk+1 делится пополам, а pk+1 пересчитывается с заменой dk+1 на dk+1/2 и так до тех пор, пока соотношение не стане меньше 10.
Функция f(x) может быть произвольного вида (не обязательно полином!)
[Б8]. Алгоритм 25б
Вызов программы
call ZR30(Fun, C, Eps1, Eps2, Eps3, Eta, Max, Cnt)
Параметры программы
Fun, Eps1, Eps2, Eps3, Eta, Max - входные параметры;
Real C[N] - входной и выходной параметр;
Real Fun(Real X) - исследуемая функция;
Real Eps1 - относительный критерий сходимости по величине промежутка, счёт прекращается, если |(pk-1 - pk-1) / pk| < Eps1;
Real Eps2 - абсолютный критерий сходимости по значениям функции;
Real Eps3, Eta - используется для кратных корней таким образом, что если |x - ci| < Eps3, то x заменяется на x+Eta;
Integer Max - максимально допустимое количество итераций;
Real С[N] - на входе: начальные приближения для корней,
на выходе: значения найденных корней;
Integer Сnt[N] - количество итераций при вычислении корня ci.
Пример
program TestZR30
use NML
implicit none
! Variables
integer, parameter:: n=6
real C(n), Eps1, Eps2, Eps3, Eta
integer Cnt(n), Max
data C/0.0, 0.0, 0.0, 0.0, 0.0, 0.0/
data Eps1/1.E-6/, Eps2/1.E-6/, Eps3/1.E-6/, Eta/0.1/, Max/100/
integer:: Count=0
integer i
! Body of TestZR30
call ZR30(Fun, C, Eps1, Eps2, Eps3, Eta, Max, Cnt)
print *, 'The result of the program:'
print 10, (i, C(i), Cnt(i), i=1,n)
print 20, Count
10 format('Root ',I1,': ',E15.7,3X,'Num.iterations =',I3)
20 format(/' Count =',I3)
contains
real function Fun(x)
real, intent(in):: x
Fun=(x-3.0)*(x-3.01)*(x-7.0)*(x-7.0)*(x+9.0)*(x+13.0)
Count=Count+1
return
end function Fun
end program TestZR30
The result of the program:
Root 1: 0.3010000E+01 Num.iterations = 16
Root 2: 0.3000000E+01 Num.iterations = 7
Root 3: 0.6999999E+01 Num.iterations = 14
Root 4: -0.9000000E+01 Num.iterations = 9
Root 5: 0.7000000E+01 Num.iterations = 5
Root 6: -0.1300000E+02 Num.iterations = 5
Count = 59
ZR40
Программа ZR40 вычисляет все корни z[1:n] полинома p(x) = a[0]*xn + ... + a[n-1]*x + a[n] степени n методом Ньютона-Миели. Полином задаётся массивом коэффициентов a[i] (i=0, 1,...,n), где a[n] - свободный член. Полином p(x) должен иметь только вещественные корни, причём все различные (не кратные). Приближения для каждого корня уточняются путём итерации, пока для двух последовательных приближений x0 и x1 не начнёт выполняться соотношение |x1-x0| ≤ eps*|x1|. При выводе корни z[i] будут упорядочены по их величинам: z[1]>z[2]>...>z[n].
[Б10]. Алгоритм 105б
Вызов программы
call ZR40(A, N, Eps, Z)
Параметры программы
A, N, Eps - входные параметры;
Integer N - степень полинома;
Real A[0:N] - коэффициенты полинома;
Real Eps - погрешность вычисления корней;
Real Z[1:N] - массив, в который помещаются найденные корни.
Пример
!Поиск вещественных корней полинома
program TestZR40
use NML
implicit none
! Variables
integer, parameter:: N=5
real A(0:N), Z(N), Eps
data A /1.0, -28.0, 74.0, 28.0, -75.0, 0.0/
data Eps /1.0E-6/
integer i
! Body of TestZR40
call ZR40(A, N, Eps, Z)
print *, 'The result of the program:'
print 20, (i, Z(i), i=1, N)
20 format('z',I1,' = ',F10.6)
end program TestZR40
The result of the program:
z1 = 25.000000
z2 = 3.000000
z3 = 1.000000
z4 = 0.000000
z5 = -1.000000
ZR44
Программа ZR44 использует корректурные формулы Брестоу и Ньютона для нахождения всех вещественных и комплексных корней полинома. Полином p(x) = a[0]*xn + ... + a[n-1]*x + a[n] задаётся массивом коэффициентов a[i] (i=0, 1,...,n), где a[n] - свободный член. В итерационном процессе коэффициенты полинома преобразуются делением их на среднее арифметическое, вычисляемое на каждом шаге. Метод итераций Брестоу-Ньютона почти всегда сходится с точностью до k-ой цифры либо к значению корня, либо к его обратному значению. Если совместные итерации по Ньютону и Брестоу не дают сходимости за it повторений, то требование к сходимости последовательно снижается на одну значащую цифру. Для каждого корня формируется параметр acc, который информирует о ходе итерационного процесса.
[Б8]. Алгоритм 30б
Вызов программы
call ZR44(A, N, it, k, U, V, Acc)
Параметры программы
A, N, it, k - входные параметры;
Integer N - степень полинома;
Real A[0:N] - коэффициенты полинома;
Integer it - максимальное количество итераций;
Integer k - желательное количество правильных значащих цифр;
Real U[1:N] - массив, содержащий вещественные части найденных корней;
Real V[1:N] - массив, содержащий мнимые части найденных корней;
Real Acc[1:N] - массив параметров, которые используются для проверки сходимости.
Пример
! Поиск вещественных и комплексных корней полинома
program TestZR44
use NML
implicit none
! Variables
integer, parameter:: N=5
real A(0:N), U(N), V(N), Acc(N)
integer it, k
data A /1.0, 1.0, -8.0, -16.0, 7.0, 15.0/
data it /100/, k /4/
integer i
! Body of TestZR44
call ZR44(A, N, it, k, U, V, Acc)
print *, 'The result of the program:'
print 20, (i, U(i), V(i), Acc(i), i=1, N)
20 format('z',I1,' = ',F10.6,' +i',F10.6,2X,'acc=',E13.6)
end program TestZR44
The result of the program:
z1 = 3.000057 +i 0.000000 acc= 0.100000E+05
z2 = -2.000029 +i -1.000046 acc= 0.100000E+05
z3 = -2.000029 +i 1.000046 acc= 0.100000E+05
z4 = -1.000000 +i 0.000000 acc= 0.100000E+05
z5 = 1.000000 +i 0.000000 acc= 0.100000E+05
!То же, для k=5
The result of the program:
z1 = 3.000000 +i 0.000000 acc= 0.100000E+06
z2 = -2.000000 +i -1.000000 acc= 0.100000E+06
z3 = -2.000000 +i 1.000000 acc= 0.100000E+06
z4 = -1.000000 +i 0.000000 acc= 0.100000E+06
z5 = 1.000000 +i 0.000000 acc= 0.100000E+06