# Faster code for solving an equation with parameters

I am new to Mathematica and wanted to solve an equation that depends on parameters. The code is below

    theta = 0.3;
v = 0.05;
alpha = 1;

q := InverseCDF[ GammaDistribution[ alpha, 1], (1 - theta)]
rk[k_, r_] := r k
rkplus1[k_, r_] := r k + 1

Lfun[k_, r_] :=
(1/rkplus1[k, r])^alpha Exp[ rk[k,r] q] CDF[ GammaDistribution[ alpha,
1/(rk[k, r] +1)],
q]

Lplustheta[k_, r_] := theta + Lfun[k, r]
logLplustheta[k_, r_] := Log[ Lplustheta[k, r]]
LHS[k_, r_] := rk[k, r] q Lfun[k, r]/Lplustheta[k, r] - logLplustheta[k, r  ]
RHS[r_] := r v
r := 1;

While[ r < 100,
{ While[ LHS[ k 0.00001, r 0.05] < RHS[r 0.05], a = k 0.00001; k++],
Print[ r 0.05, " ", v, " ", theta, " ",
Lfun[k 0.00001, r 0.05], " ", a, " ",
a/Lplustheta[ k 0.00001, r 0.05], " ",
a - (a/Lplustheta[k 0.00001, r 0.05])]; k := 1}; r++]


The code works but is extremely slow. Is there a faster way to accomplish what I trying to do?
Essentially solve RHS - LHS == 0 for k, for a range of r, and print?

• To quick things. 1) Don't use SetDelayed ( := ) to compute q, use Set instead. q is a constant and doesn't need to evaluated over and over again as you iterate through your loop. Define it once and for all with q = InverseCDF[GammaDistribution[alpha, 1], (1 - theta)]. 2) Don't write slow procedural code to solve your equation. Use one of Mathematica's built-in solvers to do it for you. – m_goldberg Aug 14 '14 at 15:11

## 1 Answer

It's your inner while loop that is causing the trouble. I refactored your code a bit also. Basically your value of k decreases after the first run. So you only have to do this whole run (over k) once. This is much faster:

θ = 0.3;
v = 0.05;
α = 1;
r = 1;
k = 1;
q = InverseCDF[GammaDistribution[α, 1], 1 - θ];
Lfun[kr_] :=
Exp[kr q]/(kr + 1)^α CDF[
GammaDistribution[α, 1/(kr + 1)], q]
LHS[kr_] :=
Block[{Z = Lfun@kr}, (kr*q*Z)/(θ + Z) - Log[θ + Z]]
While[r < 100,
If[1 == r, While[LHS[k*0.00001*r*0.05] < v*r*0.05, k++];,
While[LHS[k*0.00001*r*0.05] > v*r*0.05, --k];
++k;];
a = (k - 1)*0.00001;
σ = Lfun[k*0.00001*r*0.05];
Print[r*0.05, " ", v, " ", θ, " ", σ, " ", a, " ",
a/(θ + σ), " ", a - (a/(θ + σ))];
r++]

• Thank you very much. This works beautifully. – Shu .Anu Aug 14 '14 at 16:16