1
$\begingroup$

This question is followed up from this Question. I would like to thank Dr. Hintze and I_Mariusz for the comments and help. I am pretty new to mathematica ( I just learned it 4 days ago) so I would like to ask if there is a way to accelerate the following codes :

Clear["Global`*"]
eq = (2*S0/(y*sigma^2))^
nu*(Gamma[nu + (2*mu)/sigma^2]/Gamma[2*nu + (2*mu)/sigma^2])*
Hypergeometric1F1[nu, 2*nu + (2*mu)/sigma^2, -2*S0/(y*sigma^2)];
int = Integrate[eq, {y, K, Infinity}, GenerateConditions -> False]
int2 = int /. nu -> (-vu/2 + Sqrt[vu^2 + 8*alpha/sigma^2]/2) /. 
 vu -> (2*mu/sigma^2 - 1)

mu = 15/100;
sigma = 5/100;
S0 = 100;
K = 95;
T = 1;
F[alpha_] = int2/alpha;

A = 18.4;
n = 40;
m = 50;
S = ConstantArray[0, m + 1];
B = 1;
For[k = 0, k <= m, k++,
   S[[k + 1]] = Exp[A/2]/(2*B)*Re[F[A/(2*B)]];
   For[j = 1, j <= n + k, j++,
    S[[k + 1]] = 
S[[k + 1]] + Exp[A/2]/B*(-1)^j*Re[F[(A + 2*j*Pi*I)/(2*B)]];
]    ]
f = 0
   For[k = 0, k <= m, k++,
   f = f + Binomial[m, k]*S[[k + 1]]*2^(-m);
   ]
f

Currently, it takes me hours to run this program. Thank you so much for your time. I truly appreciate it.

$\endgroup$
8
  • 1
    $\begingroup$ I believe that if you change F[alpha_] = int2/alpha to F[alpha_] := int2/alpha this will run quite quickly. Using = in a function assignment doesn't actually assign it as a pattern (see the documentation on Set (=) versus SetDelayed (:=)). $\endgroup$ – nben Dec 15 '15 at 16:30
  • $\begingroup$ Thank you ! please let me try . $\endgroup$ – D. Nguyen Dec 15 '15 at 16:31
  • $\begingroup$ If you make that a delayed definition you need to also make int2 an explicit function of alpha. The other thing I'd suggest is to convert your For loop to a Table[ Sum[] ] structure. $\endgroup$ – george2079 Dec 15 '15 at 16:49
  • $\begingroup$ george2079 : Thanks ! , let me understand how Table is used. I've have never used it before. $\endgroup$ – D. Nguyen Dec 15 '15 at 16:51
  • 1
    $\begingroup$ The other thing you should do here is note you are frequently recomputing F for the same argument. Define F so it remembers: F[alpha_] := F[alpha] = int2[alpha]/alpha; ( or restructure your loop so you aren't repeating calculations , that inner For loop just adds one new term for each k.) $\endgroup$ – george2079 Dec 15 '15 at 16:57
3
$\begingroup$

This takes about a minute.

(* with int2 redefined as : int2[alpha_] = .. *)
F[alpha_] := F[alpha] = int2[alpha]/alpha;
terms = Monitor[ 
  Exp[A/2]/B Table[ (-1)^j*Re[F[(A + 2*j*Pi*I)/(2*B)]] ,
       {j, 0, n + m + 1}], {j}]
S = Accumulate[terms][[n ;;]]
f = Sum[Binomial[m, k]*S[[k + 1]]*2^(-m), {k, 0, m}]

I'd encourage you to study this carefully to be sure it is faithful to your calculation. (I'm sure its close but I could have shifted an index or some such doing it quickly)

$\endgroup$
0

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.