Calculating Exponent Using Taylor Series To Given Precision

[Bravus]

Homework Statement

The course is Computational Physics, but in a sense this is a pretty straight computer science or even mathematical challenge.

The first part of the assignment - the relatively easy part - was to write a Fortran program to take two variables - the number to which e should be raised, and the number of iterations of the Taylor series to use in approximating it.

The Taylor series includes the factorial of the iteration number, so a subroutine had to be written to calculate the factorial.

Got all that up and running, without too much trouble.

Homework Equations

e^x ≈ Ʃ xⁱ/i! where x is the number to evaluate. i runs from 0 to n-1, where n is the number of iterations specified.

And then there's the definition of the factorial...

The Attempt at a Solution

The second, harder (bonus marks) part of the question is to write a program that takes the same input for the number to be exponentiated, but takes another number, dx. The program should keep evaluating the Taylor Series approximation until it is less than dx away from the true value, obtained using Fortrans exp(x) intrinsic function.

The main problem I'm having is how to calculate the factorial in a loop when the upper bound keeps changing. I don't need you to write out detailed Fortran code, just some suggestions or pseudocode about how to do this calculation with these specifications.

 

[rcgldr]

The main problem I'm having is how to calculate the factorial in a loop when the upper bound keeps changing.

Can you show an attempt to calculate the factorial, even if not done in a subroutine, and an attempt to calculate xⁱ?

 

[Bravus]

Yeah, here's the code I use for the factorial, in the first part:

REAL FUNCTION rfactor1(iup1)
IMPLICIT NONE
INTEGER, INTENT(IN) :: iup1
INTEGER :: icurr1
REAL :: rans1
!
rans1=1
DO icurr1=1,iup1
rans1=rans1*icurr1
END DO
rfactor1=rans1
!
END FUNCTION rfactor1

The thing is, the variable iup1 is something the user enters, so it knows how many iterations to run. In the part I'm stuck on, the loop has to run instead until a condition is met. I've tried a DO WHILE loop, but it seems to run crazy.

In terms of 'attempts to calculate xⁱ', it's as simple in Fortran as x**1, but I assume you meant something more complex.

What I need is a loop within a loop: the inner one calculates the factorial of the counter variable of the outer loop. Sounds like it shouldn't be that hard to do, but it runs crazy, often getting into endless loops, or else delivering silly answers that are large and negative or on the order of 10^-34.

Here's my most recent attempt:

rexact2=exp(rnum2)
icount2=1 (Taylor series starts at 0 but next line takes care of that)
rapprox2=1 (Preload with zeroth term)
rdiffer2=100 (bit of a kludge, but just designing the loop to start initialised and run at least once)
DO WHILE (rdiffer2>rdx2)

DO ix=1,icount2
rxfact2=rxfact2*icount2
END DO

rapprox2=rapprox2+((rnum2**icount2)/(rxfact2))
icount2=icount2+1

END DO

rnum2 is the value of x, the number entered. rdx2 is the difference number entered - i.e. the goal is to get the estimate, rapprox2, to differ from rnum2 by less than rdx2.

A few dodgy things there mathematically, meant to get the code to run, but still not sure why it's running crazy...

 

[rcgldr]

This seems a bit more complicated that it needs to be, and I'm not sure why you chose those names for the variables. You should use double precision variables to reduce truncation error when summing up a large count of numbers.

Here is an example loop for factorial:

double precision nfac
...
      nfac = 1
      do i = 2,n
          nfac = nfac * i
      end do

Note that you could include calculation of xⁱ in this same loop as well calculation of the sum. Can you try this and show what that code would look like? You could also just calculate a summation term directly:

double precision sum, term, x
...
      x = ?
      sum = 1
      term = 1
      do i = 2,n
          term = term * ?
          sum = sum + term
      end do

 

[Bravus]

Cool, will have a bash and get back to you with results.

The variable names are because it's an assignment in a course and the lecturer/tutor wants that style... I agree, I'd far rather use simpler ones.

Thanks.

 

[Bravus]

Yeah, using double precision makes sense to avoid issues where the factorials get large, which they do pretty quickly...

 

[Bravus]

Aha! Yes, that'll work. Didn't get what you did... and then I did! Was thinking I had to run a loop from 1 to i each time to calculate the factorial, but of course I can just keep building on the factorial each time.

And yeah, forget the lecturer... I'll go for simpler variables, stop confusing myself!

 

[Bravus]

The challenge, though, is that I don't know n. There is no specific number of times to work through the loop, it just needs to run until the difference between the approximation and the actual value is smaller than a certain number.

 

[Bravus]

This works:

rexact2=exp(x)
i=1
dxfact2=1
rdiffer2=100
!
DO WHILE (rdiffer2>dx)
!
dxfact2=dxfact2*i
rapprox2=rapprox2+((x**i)/(dxfact2))
rdiffer2=abs(rapprox2-rexact2)
i=i+1
!
END DO
!

Thanks heaps, well done!

 

[rcgldr]

You can combine the xⁱ and n! calculations into one term:

      rexact2 = exp(x)
      dxterm = 1
      rapprox2 = 1
      rdiffer2 = 100
      i = 1
!
      DO WHILE (rdiffer2 > dx)
!
          dxterm = dxterm * x / i
          rapprox2 = rapprox2 + dxterm
          rdiffer2 = abs(rapprox2 - rexact2)
          i = i + 1
!
      END DO