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?
=================
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.
– N.J.Evans
Aug 11 '15 at 15:10
4
This code doesn't execute for me - something to do with the Boxes. Could you give us the InputForm of TCJS[T, k] instead?
– Patrick Stevens
Aug 11 '15 at 15:11
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.
– MarcoB
Aug 11 '15 at 15:40
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.
– J. M.♦
Aug 12 '15 at 1:14
The code seems fixed now, thanks to m_goldberg.
– Michael E2
Aug 12 '15 at 14:12
=================
1 Answer
1
=================
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} *)