0
$\begingroup$

I have written a following mathematica code which involve solution to nonlinear equations. It produce results according to my expectation but it take more than a day in execution which is very annoying. I am not expert in coding but I think my way of solving is inefficient. Can somebody help me what should I change to speedup the execution time without effecting the final results. By the way, I am using mathematica 11.3 on dell core i5, 4 GB ram machine.

1/0 // Quiet; 
α = 2.3026*(10^-3);
rho = 2.2*10^-12;
MFD = 4.2;
w0 = MFD/2;
R = List[8, 6, 4, 2, 1];
T0c = 25;
T0 = 25 + 273.15;
surrho = 1.184;  
surμ = 1.837*10^-5; 
surυ = 1.552*10^-5;  
surCp = 1.006*10^3; 
surKK = 0.02624;          
surg = 9.8;                   
 surβ = 3380*10^-6;
Array[TStack, 5];
For[x = 1, x < 6, x++,
 If[R[[x]] > MFD, Rp = MFD, Rp = R[[x]]]; 
  S0 = -(2/(Pi*(2.1*10^-6)^2))*((3.998*10^3)/(Exp[-(2*((62.5/
               w0)^2))] - 1));
 
 S = S0*Exp[-2*(r*R[[x]]/w0)^2]; 
 A0 = α*S*10^-18;
 
 surR   = R[[x]]*10^-6;        
 surκ = 
  surKK/(surCp*
     surrho);                                                         \
                                                      
 Pr = 0.71; 
 Tc = T0 + 0.000000002; 
 TcNew = Tc + 0.000000001;
 While[Abs[(((TcNew - Tc)/Tc)*100)] > 10^-10, Print["Tc"]; 
  Print[N[Tc, 25]]; Print["TcNew"]; Print[N[TcNew, 25]]; Tc = TcNew;
  Tf = (Tc + T0)/2;
  Gr = 8*surg*(surR^3)*(surυ ^-2)*(1/Tf)*(Tc - T0); 
  Print["Gr"]; Print[Gr];  
   Nu1 = (0.6 + ((0.387*((Gr*Pr)^(1/
              6)))/((1 + ((0.559/Pr)^(9/16)))^(8/27))))^2 ;
  Hc = (1/2)*Nu1*surKK*(1/surR); Print["Hc"]; Print[Hc];
  ξ = 0.85; 
  σ = 5.669*10^-8;
  Hr = (ξ*σ*((Tc^4) - (T0^4)))/(Tc - T0); 
  H = Hc + Hr; 
  H = H*10^-12 ;
  KK = (2.14749 - (298.76*Tc^(-1)) + ((20.72*10^3)*
        Tc^(-2)) - ((0.54*10^6)*Tc^-3))*10^-6;
  Cp = (1/
      60.0843)*(55.98 + ((15.40*10^(-13))*Tc) - ((14.4*10^(-5))*
        Tc^(-2)));
  κ = KK/(rho*Cp); 
  τ := κ*t;
  B0 = A0/KK ;
  h = (H/KK);
  val = (αn*R[[x]])*BesselJ[1, αn*R[[x]]] - (h*R[[x]])*
     BesselJ[0, αn*R[[x]]];
  roots = Solve[val == 0 && 0 < αn < 1, αn][[All, 1, 2]];
  Print["roots"];
  Print[roots];
  
  cn = ((2/(Rp^2))*(1/((BesselJ[
             0, (roots[[1 ;; Length[roots]]]*Rp)]^2) + (BesselJ[
             1, (roots[[1 ;; Length[roots]]]*Rp)]^2))))^0.5;  
  b1 = 2/((w0*roots[[1 ;; Length[roots]]])^2);  
  dn = Abs[
    B0*cn*((w0/
         2)^2)*(((BesselJ[0, (roots[[1 ;; Length[roots]]]*Rp)]*
           Exp[-2*((Rp/w0)^2)]) - 1) + 
       Sum[(((-1)^l)/((2^(l + 1))*(l!)*Gamma[l + 2]))*((Gamma[l + 1] -
              Gamma[l + 1, 
              b1*((roots[[1 ;; Length[roots]]]*Rp)^2)])/(2*(b1^(l + 
                 1)))), {l, 0, Infinity}])]; 
  Tn = (dn/(roots[[1 ;; Length[roots]]]^2))*(1 - 
      Exp[-(roots[[1 ;; Length[roots]]]^2)*τ]);
  u1 = T0c + (Sum[
       Tn[[i]]*cn[[i]]*BesselJ[0, (roots[[i]]*(R[[x]]*r))], {i, 1, 
        Length[roots]}]) /. t -> Infinity; Print["u1="]; Print[u1]; 
  Print["u1 at r=0"]; Print[N[u1 /. r -> 0]]; 
  Print["u1 at r=R[[x]]"]; Print[N[u1 /. r -> 1]]; 
  TcNew = N[(u1 /. r -> 1 )] + 273.15; Print["Tc="]; Print[Tc]; 
  Print["TcNew="]; Print[TcNew]];
 Print["Heat Transfer Co-efficient = "]; Print[H];
 TStack[x] = u1; Print["x"]; Print[x]; 
 Print[TStack[x]]]; Print["TStack"]; Print[TStack[1]];
Plot[{TStack[1], TStack[2], TStack[3], TStack[4], TStack[5]}, {r, 0, 
  1}, Frame -> True, PlotRange -> {{0, 1}, All}, Axes -> {True, True},
  BaseStyle -> {FontSize -> 22, FontWeight -> Plain, 
   FontFamily -> Helvetica}, 
 PlotLegends -> 
  Placed[LineLegend[(Style[#, 22, Plain, 
        FontFamily -> Helvetica] &) /@ {"R=8 [μm]", 
      "R=6 [μm]", "R=4 [μm]", "R=2 [μm]", "R=1 [μm]"},
     LegendMarkerSize -> 50, LegendLayout -> {"Row", 3}], {0.7, 0.6}],
  AxesOrigin -> {0, 25}, AxesStyle -> Dashed, 
 FrameLabel -> {"r/R (Normalized Radius)", 
   "T \!\(\*SuperscriptBox[\([\), \(o\)]\)C]"}, 
 PlotStyle -> {{Orange, Thick}, {Dashed, Red, 
    Thick}, {Dashing[{0.025, 0.01, 0.025, 0.01}], Purple, 
    Thick}, {Black, Dashing[{Large}], Thick}, {Blue, Dashing[Tiny], 
    Thick}}, ImageSize -> 700, Filling -> Axis]
```
$\endgroup$
4
  • 3
    $\begingroup$ Are all of these lines of code really necessary to produce a minimal working example of the problem you are facing? Remember that this is a Q&A site, not a free code-improving service :) $\endgroup$ Aug 21, 2020 at 9:07
  • $\begingroup$ Maybe you could also post the mathematical form of your equations in LaTeX? $\endgroup$
    – flinty
    Aug 21, 2020 at 11:05
  • $\begingroup$ @MariusLadegårdMeyer you are right but as I explained in my question the code produces result as expected but it takes very long execution time and I don't know which part of my code is inefficient. By the way I will work on it to produce a minimal working example having the same problem. $\endgroup$ Aug 21, 2020 at 13:06
  • 1
    $\begingroup$ A good start would be to find which parts/lines/commands/functions in the code take a long time to run. You can use Timing[ ] or AbsoluteTiming[ ] to help spot the slow areas. $\endgroup$
    – bill s
    Aug 21, 2020 at 13:52

1 Answer 1

3
$\begingroup$

In every line of code where u1 and TStack are using as output we need to add option /.Infinity -> 50 (we can play with it by using less then 50 or more then 50). With this option your code run few second on my machine:

α = 2.3026*(10^-3);
rho = 2.2*10^-12;
MFD = 4.2;
w0 = MFD/2;
R = List[8, 6, 4, 2, 1];
T0c = 25;
T0 = 25 + 273.15;
surrho = 1.184;
surμ = 1.837*10^-5;
surυ = 1.552*10^-5;
surCp = 1.006*10^3;
surKK = 0.02624;
surg = 9.8;
surβ = 3380*10^-6;
Array[TStack, 5];

For[x = 1, x < 6, x++, If[R[[x]] > MFD, Rp = MFD, Rp = R[[x]]];
  S0 = -(2/(Pi*(2.1*10^-6)^2))*((3.998*10^3)/(Exp[-(2*((62.5/
                w0)^2))] - 1));
  S = S0*Exp[-2*(r*R[[x]]/w0)^2];
  A0 = α*S*10^-18;
  surR = R[[x]]*10^-6;
  surκ = surKK/(surCp*surrho);
  Pr = 0.71;
  Tc = T0 + 0.000000002;
  TcNew = Tc + 0.000000001;
  While[Abs[(((TcNew - Tc)/Tc)*100)] > 10^-10, Print["Tc"];
   Print[N[Tc, 25]]; Print["TcNew"]; Print[N[TcNew, 25]]; Tc = TcNew;
   Tf = (Tc + T0)/2;
   Gr = 8*surg*(surR^3)*(surυ^-2)*(1/Tf)*(Tc - T0);
   Print["Gr"]; Print[Gr];
   Nu1 = (0.6 + ((0.387*((Gr*Pr)^(1/
               6)))/((1 + ((0.559/Pr)^(9/16)))^(8/27))))^2;
   Hc = (1/2)*Nu1*surKK*(1/surR); Print["Hc"]; Print[Hc];
   ξ = 0.85;
   σ = 5.669*10^-8;
   Hr = (ξ*σ*((Tc^4) - (T0^4)))/(Tc - T0);
   H = Hc + Hr;
   H = H*10^-12;
   KK = (2.14749 - (298.76*Tc^(-1)) + ((20.72*10^3)*
         Tc^(-2)) - ((0.54*10^6)*Tc^-3))*10^-6;
   Cp = (1/
       60.0843)*(55.98 + ((15.40*10^(-13))*Tc) - ((14.4*10^(-5))*
         Tc^(-2)));
   κ = KK/(rho*Cp);
   τ := κ*t;
   B0 = A0/KK;
   h = (H/KK);
   val = (αn*R[[x]])*BesselJ[1, αn*R[[x]]] - (h*R[[x]])*
      BesselJ[0, αn*R[[x]]];
   roots = 
    NSolve[val == 0 && 0 < αn < 1, αn][[All, 1, 2]];
   Print["roots"];
   Print[roots];
   cn = ((2/(Rp^2))*(1/((BesselJ[
              0, (roots[[1 ;; Length[roots]]]*Rp)]^2) + (BesselJ[
              1, (roots[[1 ;; Length[roots]]]*Rp)]^2))))^0.5;
   b1 = 2/((w0*roots[[1 ;; Length[roots]]])^2);
   dn = Abs[
     B0*cn*((w0/
          2)^2)*(((BesselJ[0, (roots[[1 ;; Length[roots]]]*Rp)]*
            Exp[-2*((Rp/w0)^2)]) - 1) + 
        Sum[(((-1)^l)/((2^(l + 1))*(l!)*
              Gamma[l + 2]))*((Gamma[l + 1] - 
              Gamma[l + 1, 
               b1*((roots[[1 ;; Length[roots]]]*Rp)^2)])/(2*(b1^(l + 
                  1)))), {l, 0, Infinity}])];
   Tn = (dn/(roots[[1 ;; Length[roots]]]^2))*(1 - 
       Exp[-(roots[[1 ;; Length[roots]]]^2)*τ]);
   u1 = T0c + (Sum[
        Tn[[i]]*cn[[i]]*BesselJ[0, (roots[[i]]*(R[[x]]*r))], {i, 1, 
         Length[roots]}]) /. t -> Infinity; Print["u1="]; Print[u1];
   Print["u1 at r=0"]; Print[N[u1 /. {r -> 0, Infinity -> 50}]];
   Print["u1 at r=R[[x]]"]; Print[N[u1 /. {r -> 1, Infinity -> 50}]]; 
   TcNew = N[(u1 /. {r -> 1, Infinity -> 50})] + 273.15; Print["Tc="];
    Print[Tc];
   Print["TcNew="]; Print[TcNew]]; 
  Print["Heat Transfer Co-efficient = "]; Print[H];
  TStack[x] = u1; Print["x"]; Print[x];
  Print[TStack[x]]];

Print["TStack"]; Print[TStack[1] /. Infinity -> 50];
Plot[Evaluate[{TStack[1], TStack[2], TStack[3], TStack[4], 
TStack[5]} /. Infinity -> 50], {r, 0, 1}, Frame -> True, 


PlotRange -> {{0, 1}, All}, Axes -> {True, True}, 
 BaseStyle -> {FontSize -> 22, FontWeight -> Plain, 
   FontFamily -> "Helvetica"}, 
 PlotLegends -> {"R=8 [μm]", "R=6 [μm]", "R=4 [μm]", 
   "R=2 [μm]", "R=1 [μm]"}, AxesOrigin -> {0, 25}, 
 AxesStyle -> Dashed, 
 FrameLabel -> {"r/R (Normalized Radius)", 
   "T \!\(\*SuperscriptBox[\([\), \(o\)]\)C]"}, 
 PlotStyle -> {{Orange, Thick}, {Dashed, Red, 
    Thick}, {Dashing[{0.025, 0.01, 0.025, 0.01}], Purple, 
    Thick}, {Black, Dashing[{Large}], Thick}, {Blue, Dashing[Tiny], 
    Thick}}, ImageSize -> 700, Filling -> Axis]

Figure 1

$\endgroup$
2
  • $\begingroup$ Thank You Alex. You just made my day. $\endgroup$ Aug 22, 2020 at 2:41
  • $\begingroup$ @UbaidUllah You are welcome! $\endgroup$ Aug 22, 2020 at 10:49

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.