Библиотека численных методов на языке Fortran 90. Вычисление корней функции
Автор cайта:
Владимир
Потемкин

fortran-90@yandex.ru
Язык программирования Фортран

Вычисление корней функции

>Zeroin. Вычисление корня функции f(x) в заданном интервале [a,b].

>ZR10. Вычисление корней функции f(x) при помощи изучения поведеня знака функции в интервале [a,b].

>ZR30. Вычисление корней функции f(x) методом парабол (Мюллера).

>ZR40. Метод Ньютона-Миели определения вещественных корней полинома.

>ZR44. Метод Брестоу-Ньютона определения вещественных и комплексных корней полинома.

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

 

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

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 - количество найденных корней.

 

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

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.

[Б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 - относительный критерий сходимости итераций;
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

 

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

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

 

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

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

 

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