Error in an algorithm for continued fractions
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello everyone,
I use Matlab R2020a to create an algorithm that gives all the digits of the continued fraction expansion of a real number (I know there is a function that does it more or less, but I haven't been able to use it).
function A=frac_continu(r,nmax)
% r is a real number
% nmax is the "limit rank"
A=zeros(1,nmax);
j=1;
alpha=r;
while j<= nmax && alpha-floor(alpha) ~= 0
A(j)=floor(alpha);
alpha=1/(alpha-A(j));
j=j+1;
end
end
What bothers me is that when I apply it to different real numbers, I don't get the right results from a certain rank.
For example, if I take the golden ratio (i.e. r = (1+sqrt(5))/2 and nmax=50), I should only get 1's, but the outgoing vector eventually gives me two 2's, etc.
So I wanted to know if the problem is algorithmic or if it's just due to a rounding problem or something else.
Thanks in advance!
Réponse acceptée
Rik
le 15 Fév 2021
I suspect you are indeed running into rounding issues. Using an awful approach to unwind this continued fraction, I found this for nmax=30:
r = (1+sqrt(5))/2;
nmax=30;
fprintf('%.20f (true value)\n%.20f (approximate value)\n',r,unwind(frac_continu(r,30)))
function A=frac_continu(r,nmax)
% r is a real number
% nmax is the "limit rank"
A=zeros(1,nmax);j=1;alpha=r;
while j<= nmax && alpha-floor(alpha) ~= 0
A(j)=floor(alpha); alpha=1/(alpha-A(j)); j=j+1;
end
end
function v=unwind(A)
%this is a TERRIBLE idea, don't ever use this, even a nested anonymous function is a better idea
if numel(A)==1,v=A;else,str=[sprintf('%d',A(1)),sprintf('+1/(%d',A(2:end)),repmat(')',1,numel(A)-1)];v=eval(str);end
end
After that many inversions, the rounding errors just keep on building up in your alpha, as you don't recalculate it from the A so far.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Sparse Matrices dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!