Coding a recursion in MATLAB

Question

I have been stuck trying trying to write a MATLAB algorithm that computes a recursion in reverse, or that is what it seems like to me. y_n = (1/n)−10*y_n−1 for n = 1,...,30 works in MATLAB, but because of the (*10), the round-off error makes the algorithm unstable and it is useless. Just by manipulating the recursion, y_n-1 = (1/10)(1/n - y_n) will work and the round-off errors will be reduced 10 fold at each step, potentially making this a stable algorithm.

After a couple days, I still cannot understand the logic needed to code this. Evaluating at y_n-1 is really throwing me in a loop. I was able to tackle the unstable algorithm, but I cannot think of the logic to manipulate the code to make it work. My question lies with how do you code this in MATLAB? I am truly stumped.

% Evaluate the integral yn = integral from 0 to 1 of x^n/(x+10).
% Unstable algorithm:
y(1) = log(11) - log(10);
k = 30;
for n = 1:k
    y(n+1) = (1/n) - 10*y(n);
end
n_vector = 0:k;
[n_vector;y]

By manipulating the recursion, the results will be close to true values because of the bound on the error. The current output:

0.0953101798043248 
0.0468982019567523 
0.0310179804324768 
0.0231535290085650 
0.0184647099143501 
0.0153529008564988 
0.0131376581016785 
0.0114805618403582 
0.0101943815964183 
0.00916729514692832 
0.00832704853071678 
0.00763860560192312 
0.00694727731410218 
0.00745030378205516 
-0.00307446639198020 
0.0974113305864686 
-0.911613305864686 
9.17495658805863 
-91.6940103250307 
916.992734829255 
-9169.87734829255 
91698.8211019731 
-916988.165565185 
9169881.69913012 
-91698816.9496345 
916988169.536345 
-9169881695.32499 
91698816953.2869 
-916988169532.833 
9169881695328.37 
-91698816953283.7

What is expected, with the round-off errors taken care of is the results to stay between 0and1.

I omitted just to show the values you see here. Should I include? — PattyWatty27
– PattyWatty27, Commented Sep 5, 2019 at 5:31
Well, what you need is an initial condition. If you can precompute y(30) in a similar manner (but with correct numbers) you precompute y(1), then you just reverse the order of the loop (k=30:-1:1) and just play with that — Ander Biguri
– Ander Biguri, Commented Sep 5, 2019 at 8:58
What are the values you are expecting to get? I'm not sure this is a roundoff issue at all, but just the actual result when n blows up. — Mad Physicist
– Mad Physicist, Commented Sep 5, 2019 at 10:42
Oh, that's interesting. I believe the values should stay <1. — PattyWatty27
– PattyWatty27, Commented Sep 5, 2019 at 14:12

Tirain · Accepted Answer · 2019-09-06 17:30:48Z

1

This output you are getting is correct, and as pointed out in the comments by Mad Physicist, the recursive function you have should behave this way.

If you look at the behavior of the two terms, as n gets bigger the initial subtraction will have less of an effect on the 10*y(n) term. So for large n, we can ignore 1/n.

At large n we then expect each step will increase our value by roughly a factor of 10. This is what you see in your output.

As far as writing a backward recursion. By definition you need a starting value, so you would need to assume y(30) and run the recursion backward as suggested in the comments.

answered Sep 6, 2019 at 17:30

Tirain

636 bronze badges

Sign up to request clarification or add additional context in comments.

1 Comment

PattyWatty27 Over a year ago

It shouldn't behave this way. It is a propagating error that explodes. I determined that I can only show theoretically that the absolute error is within a bound when I manipulate the equation and divide by 10 at each step, I cannot show this explicitly because I cannot compare the approximation and integral computed both by MATLAB because they both have propagating errors.

PattyWatty27 · Accepted Answer · 2019-09-07 05:38:53Z

So, I was able to answer by own question. The code needed would look like this:

% This function calculates the value of y20 with a guarantee to have an 
% absolute error less than 10^-5
% The yn1 chosen to be high enough to guarantee this is n1 = 25
% Returns the value of y(20)

function [x]= formula(k)

% RECURSION APPROXIMATION
y(k) = 0;
n = k:-1:20;
y(n-1) = (1./10)*(1./n - y(n));
x = y(20);

% FURTHER:  I needed to guarantee y20 to have <= 10^-5 magnitude error
% I determined n=25 would be my starting point, approximating y25=0 and working
% backwards to n=20 as I did above.

%  y(n-1)=1/10(1/n-yn)    “exact solution”
%  (yn-1)*=1/10(1/n-(yn)*)    “approximate solution with error”
%  y(n-1)-(y(n-1))*=1/10(1/n-yn)-1/10(1/n-(yn)*)   calculating the error
%              = 1/10((yn)*-yn)
%  So,
%  E(n-1)=1/10(En)
%  E(n-2)=1/100(E(n-1))
%  E(n-3)=1/1000(E(n-2))
%  E(n-4)=1/10000(E(n-3))
%  E(n-5)=1/100000(E(n-4)) ⇒ 10^(-5)
%  En20=(10^-5)En25
%  Therefore, if we start with n1=25, it guarantees that y20 will have 10^-5 magnitude of % the initial propagating error.

Collectives™ on Stack Overflow

Coding a recursion in MATLAB

2 Answers 2

1 Comment

Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

2 Answers 2

1 Comment

Comments

Your Answer

Sign up or log in

Post as a guest

Related