00001 program extrapolation
00002 implicit none
00003 integer :: i, k, n = 6
00004 real(8) :: xi(4),yi(4), yin(4,6), yex(4)
00005 real(8) :: v
00006 real(8), external :: lagint_4th
00007
00008 open(1,file='data.dat',status='old')
00009 do i = 1, 4
00010 read(1,*) xi(i), (yin(i,k), k = 1, n)
00011 end do
00012 close(1)
00013
00014 v = 0.0
00015
00016 do k = 1, n
00017 yi(:) = yin(:,k)
00018 yex(k) = lagint_4th(xi,yi,v)
00019 end do
00020
00021 write(6,*) (yex(k), k = 1, n)
00022 open(1,file='output_4th.dat',status='unknown')
00023 do i = 1, 4
00024 write(1,'(1p,9e18.10)')xi(i), &
00025 & (dabs(1.0d0 - yin(i,k)/yex(k))*100.0, k = 1, n)
00026 end do
00027 close(1)
00028 open(1,file='output.dat',status='unknown')
00029 do i = 1, 3
00030 write(1,'(1p,9e18.10)') &
00031 & dsqrt((xi(i)/xi(4))**2 - 1.0d0)*xi(4), &
00032 & (dabs(1.0d0 - yin(i,k)/yin(4,k))*100.0, k = 1, n)
00033 end do
00034 close(1)
00035
00036 end program extrapolation
00037
00038 function lagint_4th(x,y,v)
00039 implicit none
00040 real(8) :: lagint_4th
00041 real(8) :: x(4),y(4), v
00042 real(8) :: dx12, dx13, dx14, dx23, dx24, dx34
00043 real(8) :: dx21, dx31, dx32, dx41, dx42, dx43
00044 real(8) :: xv1, xv2, xv3, xv4, wex1, wex2, wex3, wex4
00045
00046 dx12 = x(1) - x(2)
00047 dx13 = x(1) - x(3)
00048 dx14 = x(1) - x(4)
00049 dx23 = x(2) - x(3)
00050 dx24 = x(2) - x(4)
00051 dx34 = x(3) - x(4)
00052 dx21 = - dx12
00053 dx31 = - dx13
00054 dx32 = - dx23
00055 dx41 = - dx14
00056 dx42 = - dx24
00057 dx43 = - dx34
00058 xv1 = v - x(1)
00059 xv2 = v - x(2)
00060 xv3 = v - x(3)
00061 xv4 = v - x(4)
00062 wex1 = xv2*xv3*xv4/(dx12*dx13*dx14)
00063 wex2 = xv1*xv3*xv4/(dx21*dx23*dx24)
00064 wex3 = xv1*xv2*xv4/(dx31*dx32*dx34)
00065 wex4 = xv1*xv2*xv3/(dx41*dx42*dx43)
00066
00067 lagint_4th = wex1*y(1) + wex2*y(2) + wex3*y(3) + wex4*y(4)
00068
00069 end function lagint_4th
