Comparing the value of sum of Taylor series with intrinsic function in Fortran 95

Question

I want to compare the value of sum of Taylor series with intrinsic function in Fortran 95. However, the answer for Taylor series is wrong. Below is the code I wrote

program taylor_sin
    implicit none
    
    integer, parameter :: dp = selected_real_kind(15, 307)
    integer :: i
    real(dp) ::  x_rad , term, sum, fact, y, delta
    
   write(*,*) "Enter the value of x in degree:"
    read(*,*) y
        x_rad = y * (3.14159265358979323846 / 180.0)
    
    ! Initialize the sum and the factorial
    sum = 0.0_dp
    fact = 1.0_dp
    delta = 10.0E-8
    
    ! Calculate the Taylor series expansion
    do 
     i = 1  
        term = ((-1)**i) * (x_rad**(2*i + 1)) / fact
        fact  = fact*(2*i + 1)
        sum = sum + term
       if( (sum - sin(y) ) < delta)exit
        
        
    end do
    
    
    ! Print the result
  write(*,*) " from the Taylor series"
    write(*,*) "sin(", y, ") = ", sum

        ! Print the result
    write(*,*) "intrinsic value for sin x", sin(y)
    
end program taylor_sin

This is the executable

I expect the value for both Taylor and intrinsic to be near similar

Answer 1

With sin(x), you should exploit the fact that you have an odd function and it is periodic. That is sin(-x) = -sin(x) and you can fold the argument into the range 0 <= x <= 360. As @lastchance indicates, you probably want to recursively compute the next term in the series instead of explicitly computing the factorial and the value of x raised to some integer value. Here's how I would do it with only very limited testing.

program taylor_sin

   implicit none
    
   integer, parameter :: dp =kind(1.d0)
   integer :: i, q, s
   real(dp) :: x , xs, sn, y, true
   real(dp), parameter :: deg2rad = atan(1.d0) / 45
    
   write(*,*) "Enter the value of x in degree:"
   read(*,*) y
   !
   ! Sine is an odd function.
   !
   s = 1
   if (y < 0) then
      y = -y
      s = -1
   end if
   !
   ! Sine is periodic.
   !
   y = mod(y,360._dp)
   true = sin(y * deg2rad)
   !
   ! Find quadrant. Argument reduction into 0 <= y < 90
   !
   if (y < 90) then
      s = 1    ! Do nothing!
   else if (y < 180) then
      y = y - 90
   else if (y < 270) then
      s = -s  ! Reflect through axis
      y = y - 180
   else
      s = -s  ! Reflect through axis
      y = y - 270
   end if
   !
   ! Convert from degree to radian
   !
   x = y * deg2rad
   !
   ! Use recursion to compute sine.
   ! y = sum(sn)
   ! n = 0: s0 = x
   ! n > 0: sn = - sn * x**2 / ((2 * i + 2) * (2 * i + 3))
   !
   xs = x * x / 2    ! 2 from the (2 * i + 2) term
   sn = x            ! s0 = x
   y = 0             ! Accumulator
   i = 0
   do
      y = y + sn
      sn = - sn * xs / ((i + 1) * (2 * i + 3))
      write(*,'(I2,1X,3ES23.15,ES12.4)') i, true, s * y, sn, abs(true - s * y)
      if (abs(sn) < epsilon(sn)) exit
      i = i + 1  
    end do
    
end program taylor_sin

Answer 2

Let's see. The Taylor series for sin(x) is usually written

x - x³/3! + x⁵/5! - x⁷/7! + ...

But the real problem with this is that the factorial grows very fast. Instead, write the series as ...

x - x³/(3.2.1) + x⁵/(5.4.3.2.1) - x⁷/(7.6.5.4.3.2.1) + ...

So,

first term = x

second term = -x²/(3.2) * first term

third term = -x²/(5.4) * second term

and, in general,

new term = -x²/(n*(n-1)) * previous term

where n is a variable that increases by 2 each term, reflecting the exponent of x.

So, it's most efficient to calculate each term inductively from the previous one.

We need to stop summing when this term is less than some pre-defined tolerance - not when the difference from some library function is small - because such a library function may not actually be known.

If x is large then the series would take a long time to converge. However, we do know that sin(x) is periodic with period 2.pi, so we can subtract multiples of that to reduce the range.

The code below implements this as function sin1(x).

---

However, there is another way of getting sin(x) which, ultimately, only uses the first term of the Taylor series! This is to use one of the multi-angle formulae, for example:

sin(3x)=3.sin(x)-4.sin³(x)

or

sin(x)=3.sin(x/3)-4.sin³(x/3)

Thus we can keep recursing with x/3 rather than x and this gets small very quickly. Once it is small then we have the classic small-angle formula: sin(x) = x + ...

The code below implements this as function sin2(x).

module trig
   use iso_fortran_env
   implicit none

   real(real64), parameter :: PI = 3.141592653589793238

contains

   real(real64) function sin1( arg )
      real(real64), intent(in) :: arg            ! argument in radians
   
      real(real64), parameter :: TOLERANCE = 1.0e-10
      real(real64) x, term
      integer n
   
      x = modulo( arg, 2 * PI )                  ! limit the range a bit
   
      n = 1;   term = x;   sin1 = term           ! first term
      
      do while ( abs( term ) > TOLERANCE )
         n = n + 2
         term = -term * x * x / n / ( n - 1 )    ! inductively get next term
         sin1 = sin1 + term                      ! add to series
      end do
   
   end function sin1

   !-------------------------------------------------------

   recursive real(real64) function sin2( arg ) result( res )
      real(real64), intent(in) :: arg            ! argument in radians
   
      real(real64), parameter :: TOLERANCE = 1.0e-10
      real(real64) x, S
   
      x = modulo( arg, 2 * PI )                  ! limit the range a bit
      if ( abs( x ) > TOLERANCE ) then
         S = sin2( x / 3 )
         res = S * ( 3 - 4 * S * S )             ! recursion, using sin(3x) = 3 sin(x) - 4 sin(x)^3
      else
         res = x
      end if
   
   end function sin2

end module trig

!***********************************************************************

program test
   use iso_fortran_env
   use trig
   implicit none

   integer degrees
   real(real64) radians

   do degrees = -500, 500, 10
      radians = degrees * PI / 180.0
      print "( 1x, i4, 3( 2x, f7.4 ) )", degrees, sin( radians ), sin1( radians ), sin2( radians )
   end do

end program test


 -500  -0.6428  -0.6428  -0.6428
 -490  -0.7660  -0.7660  -0.7660
 -480  -0.8660  -0.8660  -0.8660
 -470  -0.9397  -0.9397  -0.9397
 -460  -0.9848  -0.9848  -0.9848
 -450  -1.0000  -1.0000  -1.0000
 -440  -0.9848  -0.9848  -0.9848
 -430  -0.9397  -0.9397  -0.9397
 -420  -0.8660  -0.8660  -0.8660
 -410  -0.7660  -0.7660  -0.7660
 -400  -0.6428  -0.6428  -0.6428
 -390  -0.5000  -0.5000  -0.5000
 -380  -0.3420  -0.3420  -0.3420
 -370  -0.1736  -0.1736  -0.1736
 -360  -0.0000   0.0000   0.0000
 -350   0.1736   0.1736   0.1736
 -340   0.3420   0.3420   0.3420
 -330   0.5000   0.5000   0.5000
 -320   0.6428   0.6428   0.6428
 -310   0.7660   0.7660   0.7660
 -300   0.8660   0.8660   0.8660
 -290   0.9397   0.9397   0.9397
 -280   0.9848   0.9848   0.9848
 -270   1.0000   1.0000   1.0000
 -260   0.9848   0.9848   0.9848
 -250   0.9397   0.9397   0.9397
 -240   0.8660   0.8660   0.8660
 -230   0.7660   0.7660   0.7660
 -220   0.6428   0.6428   0.6428
 -210   0.5000   0.5000   0.5000
 -200   0.3420   0.3420   0.3420
 -190   0.1736   0.1736   0.1736
 -180   0.0000  -0.0000  -0.0000
 -170  -0.1736  -0.1736  -0.1736
 -160  -0.3420  -0.3420  -0.3420
 -150  -0.5000  -0.5000  -0.5000
 -140  -0.6428  -0.6428  -0.6428
 -130  -0.7660  -0.7660  -0.7660
 -120  -0.8660  -0.8660  -0.8660
 -110  -0.9397  -0.9397  -0.9397
 -100  -0.9848  -0.9848  -0.9848
  -90  -1.0000  -1.0000  -1.0000
  -80  -0.9848  -0.9848  -0.9848
  -70  -0.9397  -0.9397  -0.9397
  -60  -0.8660  -0.8660  -0.8660
  -50  -0.7660  -0.7660  -0.7660
  -40  -0.6428  -0.6428  -0.6428
  -30  -0.5000  -0.5000  -0.5000
  -20  -0.3420  -0.3420  -0.3420
  -10  -0.1736  -0.1736  -0.1736
    0   0.0000   0.0000   0.0000
   10   0.1736   0.1736   0.1736
   20   0.3420   0.3420   0.3420
   30   0.5000   0.5000   0.5000
   40   0.6428   0.6428   0.6428
   50   0.7660   0.7660   0.7660
   60   0.8660   0.8660   0.8660
   70   0.9397   0.9397   0.9397
   80   0.9848   0.9848   0.9848
   90   1.0000   1.0000   1.0000
  100   0.9848   0.9848   0.9848
  110   0.9397   0.9397   0.9397
  120   0.8660   0.8660   0.8660
  130   0.7660   0.7660   0.7660
  140   0.6428   0.6428   0.6428
  150   0.5000   0.5000   0.5000
  160   0.3420   0.3420   0.3420
  170   0.1736   0.1736   0.1736
  180  -0.0000  -0.0000  -0.0000
  190  -0.1736  -0.1736  -0.1736
  200  -0.3420  -0.3420  -0.3420
  210  -0.5000  -0.5000  -0.5000
  220  -0.6428  -0.6428  -0.6428
  230  -0.7660  -0.7660  -0.7660
  240  -0.8660  -0.8660  -0.8660
  250  -0.9397  -0.9397  -0.9397
  260  -0.9848  -0.9848  -0.9848
  270  -1.0000  -1.0000  -1.0000
  280  -0.9848  -0.9848  -0.9848
  290  -0.9397  -0.9397  -0.9397
  300  -0.8660  -0.8660  -0.8660
  310  -0.7660  -0.7660  -0.7660
  320  -0.6428  -0.6428  -0.6428
  330  -0.5000  -0.5000  -0.5000
  340  -0.3420  -0.3420  -0.3420
  350  -0.1736  -0.1736  -0.1736
  360   0.0000   0.0000   0.0000
  370   0.1736   0.1736   0.1736
  380   0.3420   0.3420   0.3420
  390   0.5000   0.5000   0.5000
  400   0.6428   0.6428   0.6428
  410   0.7660   0.7660   0.7660
  420   0.8660   0.8660   0.8660
  430   0.9397   0.9397   0.9397
  440   0.9848   0.9848   0.9848
  450   1.0000   1.0000   1.0000
  460   0.9848   0.9848   0.9848
  470   0.9397   0.9397   0.9397
  480   0.8660   0.8660   0.8660
  490   0.7660   0.7660   0.7660
  500   0.6428   0.6428   0.6428