How to make this algorithm stable

Question

I have to compute the iterative version of this formula:

f(i)=integral ( x^i/(4x+1) ) from 0 to 1

Using these formulas:

f(0)=ln(5)/4;
f(i)=1/(i+1) - 4*f(i+1);

I have tried tried the following: I calculate integral ( x^100/(4x+1) ) from 0 to 1, and store the result.Then I calculate f(i) starting from this result, using the iterative version.
But I get wrong result because the error is too big.
The error is acceptable only for i<=25.
I would like to know, why this algorithm is not stable and if there's a solution to compute the result starting from i=100 or higher.

This is the code:

function y=Integral(i,max)

if i==0
    y=1/4*log(5);
elseif i==max
    y=0.0;
else
    y=1/(i+1)-4*Integral(i+1,max);
end


end

With this function I never get an exact value because the error accumulated is too high.I get a close value (but even 3 or 4 times higher, so not acceptable) if I use i=15 and max=18.I need the stable version of this formula.

Without wanting to confirm this; is this: "f(i)=1/(i+1) - 4*f(i+1)" definition exact? Also, could you please post your code, I am a bit confused about what exactly you did and what exactly you mean by f(i) as you defined it twice... — philkark, Oct 23 '12 at 14:30
It looks more like a math homework than a programming issue.... — Oli, Oct 23 '12 at 14:35
I would agree, unless his provided code is fine and the reason for the error comes from using wrong variables and some calculations on the computer are not "exact" enough... — philkark, Oct 23 '12 at 14:36

R. Schifini · Accepted Answer · 2012-10-23T16:52:36.550

2

This recursive function should do the job without the need to store partial results on the way to 100:

function y = Integral(i)

if i==0
    y=log(5)/4;
else
    y = (-Integral(i-1) + 1/i)/4;
end
end

For recursion to work you need to start from i=100 and then call the function with i-1 until it reaches i=0.

Integral(100) will give the final answer without the need to store partial results.

edited Oct 23 '12 at 16:52

answered Oct 23 '12 at 16:31

R. Schifini

9,085
2
26
32

Yes, great this works.But there's only one more problem: if i is negative what should I do? Throw an exception or return nan? Could you show me how to handle that kind of values? – Ramy Al Zuhouri Oct 23 '12 at 18:49
1

If I'm not mistaken, for negative values of i, the integral diverges to inf, so there is no need to calculate for them, just return inf. – R. Schifini Oct 23 '12 at 19:08

How to make this algorithm stable

1 Answers1