Search code examples
matlabapproximationeulers-number

Approximate Euler's Number in script

I'm writing a script that calls a function get_fact to calculate Euler's number to 10 decimal places. My function is:

function [ nfactorial ] = get_fact( n )
%input a non-negative integer
%output is factorial of that integer
for i=0:n
    nfactorial=prod(1:n);
end
end

I now want to calculate Euler's number by having n loop from 0 to infinity and taking the sum of 1/get_fact(i).

How would I write the break in order to stop loop at the 10th decimal place? I am not allowed to use Euler's number in my script.

So far I have:

for i = 0:inf
    prod((1/ get_fact(i)))
end

but this is an infinite series.

I was thinking of using

if prod((1/ get_fact(i))) > 1E-10
break
end

but I'm not sure if this is the right way to go about the problem.

Solution

  • Firstly I would rewrite your function for finding the factorial to the simpler

    function n = factorial(n)
    n = prod(1:n);
end

    The loop in your question is unnecessary as you never use the loop variable i. I wouldn't use this function though for my solution as it can be quite slow because you have to compute redundant information at each loop iteration.

    If you still want to use a for loop you need to rewrite it to

    function f = factorial(n)
    f = 1; % 0 factorial
    for i = 1:n
        f = f * i;
    end
end

    You can use the natural logarithm and the rules of logs to determine a very accurate value of e with which you can compare against. The value of e you can check against is given by x^(1 / log(x)) where x can be any positive real number other than 1, like 2. We can see this in

                                                                    Rules of logs

    Now how do we check that we've computed a value of e to 10 decimal places of accuracy. Well given that b from above is a very accurate representation of e we can compare against it to determine when we've reached an accurate solution

    x = 2; % Any positive number other than 1
c = x^(1 / log(x));
...
if (abs(e - c) < 1e-10)
    break;
end

    In my solution e is the approximate value I've computed with the infinite sum. Note: the absolute value is taken to prevent false positives when e - c is a negative number.

    Now, an efficient method for computing the infinite sum. We can exploit how the factorial is computed to not have to compute it during each iteration hugely improving the efficiency. Firstly, we need a sum variable, e in my case, to keep track of our approximate solution. Then we need another variable to keep track of the factorial, f in my case. As 0 is a funny case we'll start off with it

    e = 0;
f = 1; % 0 factorial

e = e + 1 / f;

    and now we have the first element in our infinite sum. Next we can use the infinite sum to compute a more accurate approximate to e. The factorial can be updated during each iteration with f = f * n; leading to

    for n = 1:inf
    f = f * n; % Compute new factorial
    e = e + 1 / f; % Infinite sum
    ...
end

    Now putting that altogether produces

    x = 2; % Any positive number other than 1
c = x^(1 / log(x));

e = 0;
f = 1; % 0 factorial

e = e + 1 / f;

for n = 1:inf
    f = f * n; % Compute new factorial
    e = e + 1 / f; % Infinite sum
    if (abs(e - c) < 1e-10)
        break;
    end
end