Implement a recursive formula with internal sum

Question

I need to calculate following recursion formula.

I implemented this in MATHEMATICA as follows:

But it always gives errors for $k>0$. Can someone help me to implement this?

tn = {1.2, 2.4, 3.3};
t1 = 1.2;
sk[L_, a_, t1_, k_] := 
  sk[L, a, t1, k] = 
   a/(k + 1) Sum[(Sum[(1 - t1/tn[[j]])^i, {j, 1, L}]) sk[L, a, t1, 
       k + 1 - i], {i, 1, k + 1}];
sk[L_, a_, t1_, 0] = 1;

Answer 1

It seems the problem came from defining the function in terms of k while the mathematical formulation is defined in terms of k+1. When the two were combined, it introduced the problem of including the input term in the output, which caused an infinite loop. For example

d[2] -> With[{k = 2}, Table[d[k + 1 - i], {i, k + 1}]]

{d[2], d[1], d[0]}

would infinitely replace d[2] with the {d[2], d[1], d[0]}, whereas

d[2] -> With[{k = 2}, Table[d[k - i], {i, k}]]

{d[1], d[0]}

wouldn't output itself and thus not create an infinite loop.

Therefore changing k + 1 to k in the rhs of the definition, and leaving λ as a function for the sake of generality, gives

λreps = Thread[Array[λ, 3] -> {1.2, 2.4, 3.3}];
δ[NN_?IntegerQ, α_, 0] = 1;
δ[NN_?IntegerQ, α_, k_?IntegerQ /; k > 0] := δ[NN, α, k] = α/k Sum[
 (1 - λ[1]/λ[j])^i δ[NN, α, k - i], 
 {i, k}, {j, NN}
]

and then evaluating at k = 1 gives

δ[3, 2.3, 1]

2.3 (2 - λ[1]/λ[2] - λ[1]/λ[3])

Then for the given λ and k values

Array[δ[3, 2.3, #] &, 1 + Length@λreps, 0] /. λreps

{1, 2.61364, 4.16875, 5.23767}

Answer 2

I'd proceed as follows:

delta[0, alpha_, lambda_] = 1;

delta[k_, alpha_, lambda_] := 
 delta[k, alpha, lambda] = 
  alpha/k Sum[Sum[(1 - lambda[[1]]/lambda[[j]])^i, {j, 1, Length[lambda]}] delta[k - i, alpha, lambda], {i, 1, k}]

To test it with your values:

list = {1.2, 2.4, 3.3};
a = 2.3;

delta[1, a, list]
delta[2, a, list]
delta[3, a, list]

2.61364

4.16875

5.23767

EDIT: just a comparison on the code speed

generating random data to test:

len = 6;
vars = RandomReal[{0, 100}, len];

@That Gravity Guy code:

AbsoluteTiming[
 \[Lambda]reps = Thread[Array[\[Lambda], len] -> vars];
 \[Delta][NN_?IntegerQ, \[Alpha]_, 0] = 1;
 \[Delta][NN_?IntegerQ, \[Alpha]_, 
   k_?IntegerQ /; k > 0] := \[Delta][NN, \[Alpha], 
    k] = \[Alpha]/
     k Sum[(1 - \[Lambda][1]/\[Lambda][j])^i \[Delta][NN, \[Alpha], 
       k - i], {i, k}, {j, NN}];

 out1 = Array[\[Delta][len, 2.3, #] &, 1 + Length@\[Lambda]reps, 
    0] /. \[Lambda]reps;
 ]

0.238213

My code:

AbsoluteTiming[
 delta[0, alpha_, lambda_] = 1;
 delta[k_, alpha_, lambda_] := 
  delta[k, alpha, lambda] = 
   alpha/k Sum[
     Sum[(1 - lambda[[1]]/lambda[[j]])^i, {j, 1, 
        Length[lambda]}] delta[k - i, alpha, lambda], {i, 1, k}];
 a = 2.3;
 out2 = delta[#, a, vars] & /@ Range[0, 5];
 ]

0.00049

As you can see, the bottom code is sensibly faster, moreover it scales favourably: I doubt you can run longer lists on the first code because it scales badly with size.

If you go to a list of length 6, the first code takes 4.10815 seconds, while the bottom one takes 0.000483 on my machine.

A list of length 10 gets stuck on my machine with the first code, while it runs in 0.000602 seconds on my machine. Ofc, the outputs of the two codes are the same.