Rational Approximation for e with Matlab

Question

I've written a program that calculates the best rational approximation for e with a 1, 2, 3, 4, 5, and 6 digit denominator respectively (on Matlab). For a 5-digit denominator it takes about a minute. For a 6-digit denominator it takes 2 hours. I've tried to improve the code, but I am a beginner. Suggestions are kindly appreciated.

tic
clear
p=1;
q=1;
e=exp(1);
digits=1;
trial=1;
while numel(num2str(q))<=6
    while numel(num2str(q))==digits
        p=p+1;
        if p/q>e
            q=q+1;
        end
        A(trial)=p/q;
        B(trial,1:2)=[p,q];
        trial=trial+1;
    end
    for i=1:trial-1
        error(i)=abs(A(i)-e);
    end
    [b,r]=min(error);
    p_q(digits,1:2)=B(r,1:2);
    digits=digits+1;
end
p_q
toc

This is the output:

p_q =
          19           7
         106          39
        1264         465
       25946        9545
      271801       99990
     1084483      398959

Elapsed time is 7232.153581 seconds.

Edit: Thanks to all who answered. Dennis and Jonas were right. It was the strings and growing matrices that slowed everything down. Besides I got two wrong values for p and q. The new code is this:

tic
p_best=1;
q_best=1;
e=exp(1);
for i=1:6
    for q=10^(i-1):10^i-1
        p=round(q*e);
        if abs(p/q-e)<abs(p_best/q_best-e)
            p_best=p;
            q_best=q;
        end
    end
    disp([p_best,q_best])
end
toc

The output is:

      19           7
     193          71
    1457         536
   25946        9545
  271801       99990
 1084483      398959

Elapsed time is 0.054986 seconds.

It's faster and actually correct. I assumed I knew e like Dennis said. To Joni and hardmath: I can't use continued fractions because the teacher wanted us to use brute force, but thanks.

Answer 1

Continued fractions are an easy way to calculate good rational approximations to real numbers. You can find a disussion of an algorithm to convert a floating point number into a fraction in my blog, with code in JavaScript which should be easy to convert to matlab:

function float2rat(x) {
  var tolerance = 1.0E-6;
  var h1=1; var h2=0;
  var k1=0; var k2=1;
  var b = x;
  do {
     var a = Math.floor(b);
     var aux = h1; h1 = a*h1+h2; h2 = aux;
     aux = k1; k1 = a*k1+k2; k2 = aux;
     b = 1/(b-a);
  } while (Math.abs(x-h1/k1) > x*tolerance);
  return h1+"/"+k1;
}

If you're specifically interested in e you can take a shortcut and use the regularities in its continued fraction representation to produce rational approximations that are far more accurate than the machine epsilon. This method is described here.

Answer 2

In the evening I thought a bit about the question, and though it is a bit more radical than giving pointers on the code, I have created an example on how to approach this problem the matlab way.

Judging from your code the value of e can be considered known. With this I have created a small fragment to make over and under estimations and see which is the best one.

x=1e6:1e7-1; %To guarantee that the denominator has 6 digits 
e = exp(1);
ylow = floor(x*e); %Candidates for over estimation
yhigh = ceil(x*e); %Candidates for under estimation
lowerr = e - ylow./x ;
higherr = yhigh./x - e;
lowfind = find(lowerr == min(lowerr));
highfind = find(higherr == min(higherr));
goodfractions = [ylow(lowfind),yhigh(highfind)
                  x(lowfind),x(highfind)]
errors = [lowerr(lowfind) higherr(highfind)]   

Note that the runtime is quite fast but that the errors are or order 10^-14 which indicates that unless some measures are taken rounding errors may become relevant when you try to upscale this solution one or two digits.