1
$\begingroup$
A = 10;
a1 = 1.8;
h1 = 0.4;
b1 = 0.8;
d1 = 2900;
L1 = 0.05;
σ1 = 500;
k = 0.001;
n = 4;

c1[T_] :=
  Piecewise[{
    {3, 0 < d1*T <= 100}, 
    {2, 100 < d1*T <= 200}, 
    {1,200 < d1*T <= 300}, 
    {0.5, 300 < d1*T}, 
    {0.5, True}}];

TCJS[T_, k_] := 
  A/T + c1[T]*d1*T + h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + 
    (b1/T)*σ1*Sqrt[T + L1]*
      (PDF[NormalDistribution[0, 1], k] - 
        k*Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}]); 

EQ1[T_] := 
  (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/2)*
    (PDF[NormalDistribution[0, 1], k] - 
       k*Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}])*
    ((T + 2*L1)/(T^2*Sqrt[T + L1])) - A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

T1 = T /. Solve[EQ1[T], T, Reals][[1]];

Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]
SortBy[%, First][[1]]

I used Table because I know it is the fastest looping code in Mathematica, but it takes so much more time than matlab.

Is there any way to improve the speed?

$\endgroup$
5
  • $\begingroup$ The problem is likely the PDF's, and integrals in your definition of TCJS, because these are being re evaluated every time you evaluate an entry in the table. Using = instead of := when defining TCJS may help. But you'll need to make sure any parameters you don't explicitly list as inputs to the function are set before you define your function. $\endgroup$
    – N.J.Evans
    Aug 11 '15 at 15:10
  • 4
    $\begingroup$ This code doesn't execute for me - something to do with the Boxes. Could you give us the InputForm of TCJS[T, k] instead? $\endgroup$ Aug 11 '15 at 15:11
  • $\begingroup$ Take a look at this post on How to copy code so it looks good on this site in order to be able to provide us functional code. $\endgroup$
    – MarcoB
    Aug 11 '15 at 15:40
  • $\begingroup$ We can reopen your question as soon as you fix the issues pointed out to you, as well as provide information on what you're trying to do. $\endgroup$
    – J. M.'s torpor
    Aug 12 '15 at 1:14
  • $\begingroup$ The code seems fixed now, thanks to m_goldberg. $\endgroup$
    – Michael E2
    Aug 12 '15 at 14:12
6
$\begingroup$

Evaluate the integral once and for all (cf. cdfc):

cdfc[k_] = 
  Integrate[PDF[NormalDistribution[0, 1], y], {y, k, Infinity}];

TCJS[T_, k_] := 
  A/T + c1[T]*d1*T + 
   h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + (b1/T)*σ1*
    Sqrt[T + L1]*(PDF[NormalDistribution[0, 1], k] - k*cdfc[k]);

EQ1[T_] := (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/
       2)*(PDF[NormalDistribution[0, 1], k] - 
       k*cdfc[k])*((T + 2*L1)/(T^2*Sqrt[T + L1])) - 
    A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]; // AbsoluteTiming

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(*  {0.233491, Null}  *)

But we can extend that principal to the other functions as well, provided we block the values of T and k and use the definition of cdfc above:

Block[{T, k},

 TCJS[T_, k_] = 
  A/T + c1[T]*d1*T + 
   h1*((d1*T)/2 + k*σ1*Sqrt[T + 1]) + (b1/T)*σ1*
    Sqrt[T + L1]*(PDF[NormalDistribution[0, 1], k] - k*cdfc[k]);

 EQ1[T_] = (k*σ1*h1)/(2*Sqrt[T + L1]) - ((b1*σ1)/
       2)*(PDF[NormalDistribution[0, 1], k] - 
       k*cdfc[k])*((T + 2*L1)/(T^2*Sqrt[T + L1])) - 
    A/T^2 + (d1*h1)/2 + c1[T]*d1 == 0;

 ]

T1 = T /. Solve[EQ1[T], T, Reals][[1]];

murf = Table[{TCJS[T1, k], T1, k}, {k, 0.001, n, 0.001}]; // AbsoluteTiming

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

(*  {0.067596, Null}  *)
$\endgroup$

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.