"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)^(k1)) * (t^(alpha*(k+1))))/gamma(alpha*(k+1)+1);
> if abs(Vtemp) < 1e20
> 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 = 170e6
> 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 = 170e6, 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 roundoff 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 roundoff 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 errorprone 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
