1
$\begingroup$

I need to calculate following recursion formula. enter image description here

I implemented this in MATHEMATICA as follows: enter image description here

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;
$\endgroup$
  • 1
    $\begingroup$ Can you please post your code so we can copy and paste easily? $\endgroup$ – A_user_with_NoName Oct 3 at 18:45
1
$\begingroup$

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}

$\endgroup$
  • $\begingroup$ I am confused with two sets of different outputs from here and @Fraccalo. Any idea? $\endgroup$ – Frey Oct 3 at 23:39
  • $\begingroup$ Good explanation and helpful answer! $\endgroup$ – Frey Oct 4 at 0:01
2
$\begingroup$

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.

$\endgroup$
  • 1
    $\begingroup$ I think you made a typo and forgot the i exponent in your summand. $\endgroup$ – That Gravity Guy Oct 3 at 23:48
  • $\begingroup$ I even didn't notice it @That Gravity Guy. Nice simple code and useful. $\endgroup$ – Frey Oct 4 at 0:03
  • $\begingroup$ absolutely @ThatGravityGuy, fixed it now. Thanks $\endgroup$ – Fraccalo Oct 4 at 7:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.