sum of alternating series

Hello there!

I am trying to find the sum of a sign alternating series (with a finite number of terms). The series is as follows:

\rho_K(x,A) = 1 + \sum_{n=1}^K C_n^K \exp\{(b/\lambda) (1-\lambda)^n x\},

where

C_n^K = C_{n-1}^{K-1} / (1 - (1 - lambda)^{n-1}) for n = 2,...,K
and C_1^K = - exp\{-(b/\lambda) * A\} (1 + \sum_{n=2}^K C_n^K \exp\{(b/\lambda) (1-\lambda)^{n-1} A\})
and C_1^1=-exp\{-(b/\lambda) * A\}.

I've implemented all this in MATLAB (see script below), it works, but ... The problem is that the coefficients C_1^K, C_2^K, ..., C_K^K are alternating in sign, and when lambda is small (as in 0.01 or less), and K is large (as in 100 or so) the sum explodes and returns nonsense. The sum is supposed to be between 0 and 1 (it's a probability of a certain event), no matter what K is, and it's supposed to decrease as K gets larger. The particular set of parameters I've tried is lambda = 0.01, A = 1.3, x = 0, b = 1, K = 100. I would appreciate any help to make the series more stable and easier to add. I've already tried splitting the series into two parts: one with negative coefficients, and one with positive coefficients, and treating the two separately - didn't work.

Thanks a lot in advance,

function rho_k = fn_rho_k(lambda, A, x, b, K)

% powers of (1 - lambda) from 1 thru K inclusive
lambda_arr = cumprod(ones(1, K) .* (1 - lambda));
exp_arr = exp((b/lambda) .* lambda_arr .* A);

% array of constants
C = zeros(1, K, 'double');

% start with the constant for rho1
C(1) = -1 / exp((b / lambda) * A);

for i = 2:K
    C(2:i) = C(1:(i-1)) ./ (1 - lambda_arr(1:(i-1)));
    C(1) = -(1 + dot(C(2:i), exp_arr(1:(i-1)))) / exp((b / lambda) * A);
end

rho_k = ones(size(x));

% rho_k = rho_k + dot(C, exp((b / lambda) .* lambda_arr(:) .* x));

for i = 1:K
    rho_k = rho_k + C(i) * exp((b / lambda) .* lambda_arr(i) .* x);
end

end

[Jan]

Do you have any evidence, that you have implemented the formula correctly?

Answer 1

[Jan]

When I run you program with the provide test data

lambda = 0.01, A = 1.3, x = 0, b = 1, K = 100

I get the result:

rho_k = 0.999963563929568

If you expect the result in the interval [0, 1], this looks fine.

I've stored the elements of the final sum ("for i=1:K") in an array and calculated the sum using the more precise XSum to get exactly the same result.

Please provide any example values, which reproduce your "exploding" values.

[Alex Poe]

May I ask which version of MATLAB you use? I tried running the script on R2012a and on R2012b with the same parameters, and obtained two different results (and different from yours). Try K=200, or increase it even further 300 or 400. Does it still work? For instance, this is what I get on R2012b for K=200:

>> fn_rho_k(0.01,1.3,0,1,200)

ans =

            82946736.8513228

Thanks much for your help.

[Alex Poe]

Here's another example: lambda = 0.01, A = 1.14047748184992, x = 1, b = 1, K = 100

and the result is 1.01634381235788e+22

This is using MATLAB R2012b, 64bit.

[Jan]

I get the same results with R6.5 (remove the not supported 'double' from the zeros() command), R2009a/32 and 64 and R2011b/64 und WindowsXP/32 and 7/64.

I'll try your new inputs in the evening. But it is strange already, that you get different results for 2012a and 2012b. Can you provide any correct result?

Do you have any shadowed built-in functions in your path, e.g. an M-file called exp.m? Try to remove all user-defined folders from the path and run your program again.