Thread Subject:

rounding issue

Hi!

I have to calculate a voltage raise in a constant phase element, using following code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I0 = -0.001; % pulse current
k = 0;
temp=0;

while(k<1000)
    Vtemp = (((-1)^k) * ((R*C)^(-k-1)) * (t^(alpha*(k+1))))/gamma(alpha*(k+1)+1);
    if abs(Vtemp) < 1e-20
        break
    else
        k = k+1;
        temp = temp + Vtemp;
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

This calculates the voltage depletion as a function of the time, and of R, C and alpha, which are parameters of the CPE. k typically rises up to 200.

Now, my problem is, using:
R = 470
C = 170e-6
alpha = 0.62

the function converges for time values inferior to 0.4 [s], but then gives incoherent results. I suspect this is a rounding error, due to the division of two very large numbers.

Does anyone have an idea how to circumvent this??

thanks a lot!

"G. B." <gavrilo.dot.bozovic@gmail.dot.ch> wrote in message <gq04jp$8ee$1@fred.mathworks.com>...
> Hi!
>
> I have to calculate a voltage raise in a constant phase element, using following code:
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> I0 = -0.001; % pulse current
> k = 0;
> temp=0;
>
> while(k<1000)
> Vtemp = (((-1)^k) * ((R*C)^(-k-1)) * (t^(alpha*(k+1))))/gamma(alpha*(k+1)+1);
> if abs(Vtemp) < 1e-20
> break
> else
> k = k+1;
> temp = temp + Vtemp;
> end
> end
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> This calculates the voltage depletion as a function of the time, and of R, C and alpha, which are parameters of the CPE. k typically rises up to 200.
>
> Now, my problem is, using:
> R = 470
> C = 170e-6
> alpha = 0.62
>
> the function converges for time values inferior to 0.4 [s], but then gives incoherent results. I suspect this is a rounding error, due to the division of two very large numbers.
>
> Does anyone have an idea how to circumvent this??
>
> thanks a lot!

  I'm afraid you won't like what I have to tell you, G. B. I tried out your code for the R = 470, C = 170e-6, alpha = .62, and t = 0.4 combination, and it is rather clear what the difficulty is, but it isn't evident what to do to improve things.

  You are correct that there are some large quantities being divided by other large quantities, but for these numbers that isn't the cause of the trouble. The problem lies in the reversing signs in the series being computed and the magnitudes of the Vtemp terms being subtracted from one another. At around the range from k = 20 to k = 40, the magnitude of Vtemp climbs up precipitously to somewhere around 10^9 or thereabouts, so you have very large numbers being added and subtracted sequentially to form a small final answer near one. This is a procedure guaranteed to produce large round off errors, since the round-off errors are approximately 1/10^16 times the magnitude of the terms being combined.

  For larger values of t this becomes far more extreme. For example, with t = 0.7 the Vtemp magnitude climbs to around 10^16 even though the final answer is presumably still somewhere in the neighborhood of one, and this means that there is an essentially total loss of accuracy, with the final answer being of about the same size as a number of consecutive round-off errors - the final answer would be totally obscured by the errors.

  Offhand I would say that the very nature of the series you are using is doomed to fail in such circumstances, and there is a need for a different mathematical method for computing your desired "voltage raise" or whatever it is you are calculating.

  Your function is of the form

 f(x,y) = K * SUM, k = 1 to infinity, (-x)^k/gamma(y*k)

with x > 0, and it is an inherently error-prone formula for the wrong combinations of x and y values. The values x and y here correspond to the values t^alpha/(R*C) and alpha, respectively. They are about 10 and 0.62 in the t = .7 case, and around 7 and .62 in the t = .4 case.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gq1pi9$87i$1@fred.mathworks.com>...
> ......
> Offhand I would say that the very nature of the series you are using is doomed to fail in such circumstances, and there is a need for a different mathematical method for computing your desired "voltage raise" or whatever it is you are calculating.
> ......

  I should have pointed out that your series is quite analogous to calculating exp(-x) using the Taylor expansion where -x is a large negative number such as -20:

 exp(-x) = 1/0! - x/1! + x^2/2! - x^3/3! + x^4/4! - x^5/5! + ...

The correct sum of the series is exp(-20) = 2.06e-9 but the individual terms in the series reach a maximum of (-20)^20/20! = 4.31e+7. The round-off error in computing that one term alone can be larger than the correct answer, so there is no hope of calculating a decent exp(-20) in this manner.

Roger Stafford

"G. B." <gavrilo.dot.bozovic@gmail.dot.ch> wrote in message <gq04jp$8ee$1@fred.mathworks.com>...
> Hi!
>
> I have to calculate a voltage raise in a constant phase element, using following code:
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
>
> I0 = -0.001; % pulse current
> k = 0;
> temp=0;
>
> while(k<1000)
> Vtemp = (((-1)^k) * ((R*C)^(-k-1)) * (t^(alpha*(k+1))))/gamma(alpha*(k+1)+1);
> if abs(Vtemp) < 1e-20
> break
> else
> k = k+1;
> temp = temp + Vtemp;
> end
> end
>

Was the above sum approximates an integral? How alpha is chosen? I guess it is related to a timestep? Have you analyzed what kind of error induced for the alpha you are using?

Sometime it's better to go back to the blackboard and take a new look on the original problem instead of persisting on the bad discretization.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gq2ak0$gtp$1@fred.mathworks.com>...
> "G. B." <gavrilo.dot.bozovic@gmail.dot.ch> wrote in message <gq04jp$8ee$1@fred.mathworks.com>...
> > Hi!
> >
> > I have to calculate a voltage raise in a constant phase element, using following code:
> >
> > %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> >
> > I0 = -0.001; % pulse current
> > k = 0;
> > temp=0;
> >
> > while(k<1000)
> > Vtemp = (((-1)^k) * ((R*C)^(-k-1)) * (t^(alpha*(k+1))))/gamma(alpha*(k+1)+1);
> > if abs(Vtemp) < 1e-20
> > break
> > else
> > k = k+1;
> > temp = temp + Vtemp;
> > end
> > end
> >
>
> Was the above sum approximates an integral? How alpha is chosen? I guess it is related to a timestep? Have you analyzed what kind of error induced for the alpha you are using?
>
> Sometime it's better to go back to the blackboard and take a new look on the original problem instead of persisting on the bad discretization.
>
> Bruno

I'd very much agree with Bruno. One solution, IMHO
the wrong one, is to immediately look for higher
precision. The problem is, this just pushes the problem
out a little ways. And then you need to go to higher
precision yet. And high precision arithmetic is slow.

It is better to look for a different approach. Many
times it is something as simple as scaling the problem
so that the parameter lies in the domain where the
series is nicely convergent.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gq2efk$src$1@fred.mathworks.com>...
> "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gq2ak0$gtp$1@fred.mathworks.com>...
> > Was the above sum approximates an integral? How alpha is chosen? I guess it is related to a timestep? Have you analyzed what kind of error induced for the alpha you are using?
> >
> > Sometime it's better to go back to the blackboard and take a new look on the original problem instead of persisting on the bad discretization.
> >
> > Bruno
>
> I'd very much agree with Bruno. One solution, IMHO
> the wrong one, is to immediately look for higher
> precision. The problem is, this just pushes the problem
> out a little ways. And then you need to go to higher
> precision yet. And high precision arithmetic is slow.
>
> It is better to look for a different approach. Many
> times it is something as simple as scaling the problem
> so that the parameter lies in the domain where the
> series is nicely convergent.
>
> John

  Bruno and John, what you say may be true. However, it remains a fact that the function defined by Gavrilo's infinite series is a very respectable entire analytic function, f(z), (analytic over the entire complex plane) using z = t^alpha/(R*C). For alpha = 1, it is the function

 f(z) = 1-exp(-z)

and for alpha = 2, it is

 f(z) = 1-cos(z^(1/2))

both of them entire analytic functions. The only troublesome thing about them is that their Taylor series is expanded about the wrong point. About a nearer point there should be an accurate way to compute them. The problem seems to be one of applying the wrong numerical technique for its evaluation, not something intrinsic in the function itself.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gq3bjh$jj$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <gq2efk$src$1@fred.mathworks.com>...
> > "Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <gq2ak0$gtp$1@fred.mathworks.com>...
> > > Was the above sum approximates an integral? How alpha is chosen? I guess it is related to a timestep? Have you analyzed what kind of error induced for the alpha you are using?
> > >
> > > Sometime it's better to go back to the blackboard and take a new look on the original problem instead of persisting on the bad discretization.
> > >
> > > Bruno
> >
> > I'd very much agree with Bruno. One solution, IMHO
> > the wrong one, is to immediately look for higher
> > precision. The problem is, this just pushes the problem
> > out a little ways. And then you need to go to higher
> > precision yet. And high precision arithmetic is slow.
> >
> > It is better to look for a different approach. Many
> > times it is something as simple as scaling the problem
> > so that the parameter lies in the domain where the
> > series is nicely convergent.
> >
> > John
>
> Bruno and John, what you say may be true. However, it remains a fact that the function defined by Gavrilo's infinite series is a very respectable entire analytic function, f(z), (analytic over the entire complex plane) using z = t^alpha/(R*C). For alpha = 1, it is the function
>
> f(z) = 1-exp(-z)
>
> and for alpha = 2, it is
>
> f(z) = 1-cos(z^(1/2))
>
> both of them entire analytic functions. The only troublesome thing about them is that their Taylor series is expanded about the wrong point. About a nearer point there should be an accurate way to compute them. The problem seems to be one of applying the wrong numerical technique for its evaluation, not something intrinsic in the function itself.
>
> Roger Stafford

Roger, John, Bruno,

Thanks a lot for your very professional inputs. I understand better the source of the problem.

Just in case you are interested: the time series described is one solution in the time domain for voltages in batteries containing a constant phase element, whose formula is:

Z(omega) = 1/((j*omega)^alpha * C)

I'll try to solve this in frequency domain, which shouldn't bring these errors!

Thanks again, your help is invaluable!