Taylor Series Expansion for exp(x)/sin(x) in Fortran
-
27-05-2021 - |
Вопрос
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
Решение
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.
Другие советы
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.