math - Constructing Romberg integration table -
i trying write fortran program generate romberg integration table. there's algorithm in book numerical analysis r.l.burden , j.d.faires 9th ed. in chapter 4.5. far have written this
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) i=2,n sum1=0. k=1,2**(i-2) sum1=sum1+f(a+(k-.5)*h) enddo r(2,1)=.5*(r(1,1)+h*sum1) 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. 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 this program gives following output:
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 but supposed give this:
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 the above 1 done program written in maple. program in maple is
> 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))); row 1 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 col 1 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) row 0 n do; col 0 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) so fortran program same result found maple program?
there number of differences between maple program , fortran source.
the result array of maple program dimensioned 0 n, while fortran program runs 1 n.
the fortran source never defines (calculates value for) r(3:,:) on account of fixed column indices.
given differences, shouldn't surprising results differ.
a naive, relatively direct, translation of maple source f2008 gives same result, after accounting usual vagaries of floating point arithmetic.
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) :: 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 - r(0,0) = h / 2 * (f(a) + f(b)) 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))]) col = 1, row r(row, col) = (4**col * r(row, col-1) - r(row-1, col-1)) & / (4**col - 1) end end end function romberg subroutine print_table(unit, r) integer, intent(in) :: unit real(rk), intent(in) :: r(0:,0:) integer :: row row = 0, ubound(r,1) write (unit, "(*(f13.10,1x))") r(row, :row) end 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
Comments
Post a Comment