Построение интегральной таблицы Ромберга

Я пытаюсь написать программу на Фортране для генерации таблицы интеграции Ромберга. В книге " Численный анализ" Р.Л.Бурдена и JDFaires, 9-е изд., Есть алгоритм . в главе 4.5. Пока я написал это

implicit none
integer,parameter::n=4
real::a,b,f,r(n,n),h,sum1
integer::i,k,j,m,l
open(1,file='out.txt')
a=0.
b=1.
h=b-a
r(1,1)=.5*h*(f(a)+f(b))
write(1,*)r(1,1)
do i=2,n
    sum1=0.
    do k=1,2**(i-2)
        sum1=sum1+f(a+(k-.5)*h)
    enddo
    r(2,1)=.5*(r(1,1)+h*sum1)
    do j=2,i
        r(2,j)=r(2,j-1)+(r(2,j-1)-r(1,j-1))/(4**(j-1)-1)
        write(1,*)((r(m,l),m=2,2),l=1,i)
    enddo
    h=h/2.
    do j=1,i
        r(1,j)=r(2,j)
    enddo
enddo

end

real function f(x)
implicit none
real,intent(in)::x

f=1/(1+x**2)

end function

Эта программа дает следующий вывод:

  0.750000000    
  0.774999976      0.783333302    
  0.782794118      0.785392165       3.56011134E-22
  0.782794118      0.785392165      0.785529435    
  0.784747124      0.785398126      0.785529435       7.30006976E+28
  0.784747124      0.785398126      0.785398543       7.30006976E+28
  0.784747124      0.785398126      0.785398543      0.785396457

Но это должно дать:

0.7500000000 
0.7750000000 0.7833333333 
0.7827941176 0.7853921567 0.7855294120 
0.7847471236 0.7853981253 0.7853985227 0.7853964451 
0.7852354030 0.7853981627 0.7853981647 0.7853981590 0.7853981659

Вышеуказанное сделано программой, написанной на Maple. Программа в Maple есть

>   romberg := proc(f::algebraic, a, b, N,print_table) 
 local R,h,k,row,col; 
 R := array(0..N,0..N); 

 # Compute column 0, Trapezoid Rule approximations of 
 #                   1,2,4,8,..2^N subintervals 
 h := evalf(b - a); 
 R[0,0] := evalf(h/2 * (f(a)+f(b))); 
 for row from 1 to N do; 
   h := h/2; 
   R[row,0] := evalf(0.5*R[row-1,0] + 
                     sum(h*f(a+(2*k-1)*h),k=1..2^(row-1))); 
   # Compute [row,1]:[row,row], via Richardson extrapolation 
   for col from 1 to row do; 
     R[row,col] := ((4^col)*R[row,col-1] - R[row-1,col-1]) / 
                       (4^col - 1); 
   end do; 
 end do; 

 # Display results if requested 
 if (print_table) then 
       for row from 0 to N do; 
     for col from 0 to row do; 
       printf("%12.10f ",R[row,col]); 
     end do; 
     printf("\n"); 
   end do; 
 end if; 

 return(R[N,N]); 

end proc:
f:=x->1/(1+x^2);
val:=romberg(f,0,1,4,true)

Так что же теперь делать с программой на Фортране, чтобы получить тот же результат, что и в программе Maple?

Источник

2 ответа

Есть много различий между программой клена и источником Фортрана.

  • Массив результатов программы maple имеет размер от 0 до n, а программа Fortran работает от 1 до n.

  • Источник на Фортране никогда не определяет (вычисляет значение) r(3:,:) из-за фиксированных индексов столбцов.

Учитывая эти различия, неудивительно, что результаты отличаются.

Наивный, относительно прямой перевод источника клена на F2008 дает тот же результат после учета обычных капризов арифметики с плавающей запятой.

module romberg_module
  implicit none

  integer, parameter :: rk = kind(1.0d0)

  abstract interface
    function f_interface(x)
      import :: rk
      implicit none
      real(rk), intent(in) :: x
      real(rk) :: f_interface
    end function f_interface
  end interface
contains
  function romberg(f, a, b, n) result(r)
    procedure(f_interface) :: f
    real(rk), intent(in) :: a
    real(rk), intent(in) :: b
    integer, intent(in) :: n
    real(rk) :: r(0:n,0:n)    ! function result.

    real(rk) :: h
    integer :: row
    integer :: col
    integer :: k

    h = b - a
    r(0,0) = h / 2 * (f(a) + f(b))
    do row = 1, n
      h = h / 2
      r(row, 0) = 0.5_rk * r(row-1, 0)  &
          + sum(h * [(f(a + (2 * k - 1) * h), k = 1, 2**(row-1))])
      do col = 1, row
        r(row, col) = (4**col * r(row, col-1) - r(row-1, col-1))  &
            / (4**col - 1)
      end do
    end do
  end function romberg

  subroutine print_table(unit, r)
    integer, intent(in) :: unit
    real(rk), intent(in) :: r(0:,0:)
    integer :: row
    do row = 0, ubound(r,1)
      write (unit, "(*(F13.10,1X))") r(row, :row)
    end do
  end subroutine print_table
end module romberg_module

program print_romberg_table
  use, intrinsic :: iso_fortran_env, only: output_unit
  use romberg_module
  implicit none
  real(rk), allocatable :: r(:,:)
  r = romberg(f, 0.0_rk, 1.0_rk, 4)
  call print_table(output_unit, r)
contains
  function f(x)
    real(rk), intent(in) :: x
    real(rk) :: f
    f = 1.0_rk / (1.0_rk + x**2)
  end function f
end program print_romberg_table
Источник

Ваша программа слишком длинна, чтобы ее можно было отлаживать, просто читая код (и уже более 5 лет я не пишу ни одной строки из f95), но из вашей таблицы результатов видно, что в F90 реальные переменные представлены в виде одинарной точности (32 бита) IEEE 754 числа с плавающей запятой, в то время как клен использует более высокую точность.

Фактически для 32-битной точности 0,775 становится (приблизительно) 0,77499998, а для 64-битной - 0,77500000000000002. Абсолютно неправильные цифры (например, 3.56011134E-22) могут быть вызваны недостаточным / переполнением или отменой.

Вы можете добиться 64-битной точности с помощью подходящих опций компилятора или путем указания соответствующего реального вида, см. http://fortranwiki.org/fortran/show/Real+precision

integer, parameter :: dp = kind(1.d0)
real(kind=dp) :: a
Источник
Другие вопросы по тегам