Taylor Series Expansion for exp(x)/sin(x) in Fortran

Question

I tried to write a Taylor series expansion for exp(x)/sin(x) using fortran, but when I tested my implementatin for small numbers(N=3 and X=1.0) and add them manually, the results are not matching what I expect. On by hand I calculated 4.444.., and with the program I found 7.54113. Could you please check my code and tell me if i got anything wrong.

Here is the expansion formula for e^x/sin(x) in wolframalpha: http://www.wolframalpha.com/input/?i=e%5Ex%2Fsin%28x%29

 PROGRAM Taylor
 IMPLICIT NONE


 INTEGER ::Count1,Count2,N=3
 REAL:: X=1.0,Sum=0.0
 COMPLEX ::i=(0.0,0.1)
 INTEGER:: FACT

  DO Count1=1,N,1
    DO Count2=0,N,1
       Sum=Sum+EXP(i*X*(-1+2*Count1))*(X**Count2)/FACT(Count2)
    END DO

 END DO

 PRINT*,Sum

 END PROGRAM Taylor  

 INTEGER FUNCTION FACT(n)
   IMPLICIT NONE
   INTEGER, INTENT(IN) :: n
   INTEGER :: i, Ans
   Ans = 1
   DO i = 1, n
     Ans = Ans * i
   END DO
   FACT = Ans
END FUNCTION FACT

Answer 1

The Wolfram article has the expansion formula stated using q = e**(ix) so there is a complex term. Therefore "sum" should be declared complex.

As already stated, the factorial function is simplistic. Be careful of overflow.

It is best to place your procedures into a module and "use" that module from the main program. Use as many compiler debugging options as possible. For example, gfortran, when the appropriate warning option is used, warns about the type of "sum": "Warning: Possible change of value in conversion from COMPLEX(4) to REAL(4)". If you are using gfortran try: -O2 -fimplicit-none -Wall -Wline-truncation -Wcharacter-truncation -Wsurprising -Waliasing -Wimplicit-interface -Wunused-parameter -fwhole-file -fcheck=all -std=f2008 -pedantic -fbacktrace

Since you can do this problem by hand, try outputting each step with a write statement and comparing to your hand calculation. You will probably quickly see where the calculation diverges. If it isn't clear why the calculation is different, break it down into pieces.

Answer 2

I don't see any complex terms in that expansion by Wolfram, so I'd wonder why you think you need the complex number in the exponential term. And you can't get that 1/x term the way you've programmed it. You need an x**(-1.0) term somewhere.

Your factorial implementation is rather naive, too.

I'd recommend that you forget about the loop and factorials and start with a polynomial, coefficients, and Horner's method for evaluation. Get that working and then see if you can sort out the loop.