2
$\begingroup$

I am trying to use the code developed by @Alex Trounev here: Numerical solution of an iterative equation for different kind of q. I am not sure how to incorporate two different parameters in a Do loop to achieve this.

Here's my manual try to plot Tf and cp for different q:

These are the parameters in common for any q:

ah = 342496; (*J/mol*)
R = 8.314; (*J/mol.K*)
A = 7.6*10^-38; (*s*)
b = 0.67;
x = 0.49;
T0 = 500;(*K*)
Tfinal = 350;(*K*)
Tf[0] = T0;(*Initial condition for Tf*)

Here's to implement a q=-0.1/60 (q1) and plot it:

(*q1*)
q = -0.1/60;
dt = Abs[(T0 - Tfinal)/(300*q)];(*s. This ensures n to be 300*)
n = IntegerPart[Abs[(T0 - Tfinal)/dt/q]]; (*number of steps*)
dT = dt*q; (*K*)
T[k_] := T0 + \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 1\), \(k\)]dT\);

Do[
 tau[k] = A*Exp[(x ah)/(R T[k]) + (1 - x) ah/(R Tf[k - 1])];
 
 Tf[k] = T0 + \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 
      1\), \(k\)]\(dT*\((1 - Exp[\((\(-\(
\*SubsuperscriptBox[\(\[Sum]\), \(j = 
             i\), \(k\)]\((\((dT)\)/\((q*tau[j])\))\)^b\)\))\)])\)\)\);
 
 , {k, n} (*Do[exp,{i,imax}]. k starts at 1.*)]

(*Visualization for q1*)

Tfc = Table[{T[k], Tf[k]}, {k, 
   n}]; (*Creates a table of Tk and Tfk such \
{{Tk1,Tfk1},{Tk2,Tfk2}...} from k to n*)
T[0] = T0;
cp = Table[{T[k], (Tf[k] - Tf[k - 1])/(T[k] - T[k - 1])}, {k, n}];

ListPlot[Tfc, PlotRange -> {{350, 500}, {400, 500}}, 
 ImageSize -> Medium, AxesLabel -> {"T", "Tf"}]
ListPlot[cp, PlotRange -> {{350, 500}, {-0.5, 3}}, 
 ImageSize -> Medium, AxesLabel -> {"T", "cp"}]

Here's to implement a q=-1/60 (q2) and plot it (the same except for changing q):

(*q2*)
q = -1/60;
dt = Abs[(T0 - Tfinal)/(300*q)];(*s. This ensures n to be 300*)
n = IntegerPart[Abs[(T0 - Tfinal)/dt/q]]; (*number of steps*)
dT = dt*q; (*K*)
T[k_] := T0 + \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 1\), \(k\)]dT\);

Do[
 tau[k] = A*Exp[(x ah)/(R T[k]) + (1 - x) ah/(R Tf[k - 1])];
 
 Tf[k] = T0 + \!\(
\*SubsuperscriptBox[\(\[Sum]\), \(i = 
      1\), \(k\)]\(dT*\((1 - Exp[\((\(-\(
\*SubsuperscriptBox[\(\[Sum]\), \(j = 
             i\), \(k\)]\((\((dT)\)/\((q*tau[j])\))\)^b\)\))\)])\)\)\);
 
 , {k, n} (*Do[exp,{i,imax}]. k starts at 1.*)]


(*Visualization for q2*)

Tfc = Table[{T[k], Tf[k]}, {k, 
   n}]; (*Creates a table of Tk and Tfk such \
{{Tk1,Tfk1},{Tk2,Tfk2}...} from k to n*)
T[0] = T0;
cp = Table[{T[k], (Tf[k] - Tf[k - 1])/(T[k] - T[k - 1])}, {k, n}];

ListPlot[Tfc, PlotRange -> {{350, 500}, {400, 500}}, 
 ImageSize -> Medium, AxesLabel -> {"T", "Tf"}]
ListPlot[cp, PlotRange -> {{350, 500}, {-0.5, 3}}, 
 ImageSize -> Medium, AxesLabel -> {"T", "cp"}]

Question:

What I want is to be able to do each plot for different q's in a single plot in a automatic manner without having to repeat the code at each q. The q values I want to plot are q1=-0.1/60,q2=-1/60,q1=-10/60, q1=-80/60. How can I do this?

$\endgroup$
1
$\begingroup$
  1. Wrap everything in a Module (added bonus: localize variables and definitions)
  2. Use it to define a function that takes q as a parameter
  3. Build a Table that calls your function with the appropriate values of q
  4. Done!
ClearAll[function]
function[q_, T0_: 500, Tfinal_: 350, ah_: 342496] :=
 Module[{R, A, b, x, dt, n, dT, T, Tf, tau, Tfc, cp},
  R = 8.314;(*J/mol.K*)
  A = 7.6*10^-38;(*s*)
  b = 0.67;
  x = 0.49;
  Tf[0] = T0;(*Initial condition for Tf*)
  
  dt = Abs[(T0 - Tfinal)/(300*q)];
  n = IntegerPart[Abs[(T0 - Tfinal)/dt/q]];
  dT = dt*q;
  T[k_] := T0 + Sum[dT, {i, 1, k}];
  
  Do[
   tau[k] = A Exp[(x ah)/(R T[k]) + (1 - x) (ah/(R Tf[k - 1]))];   
   Tf[k] = T0 + Sum[dT (1 - Exp[-Sum[(dT/(q*tau[j]))^b, {j, i, k}]]), {i, 1, k}],
   {k, n}
  ];
  
  Tfc = Table[{T[k], Tf[k]}, {k, n}]; T[0] = T0; 
  cp = Table[{T[k], (Tf[k] - Tf[k - 1])/(T[k] - T[k - 1])}, {k, n}];
  {Tfc, cp}
 ]

Then you can call function with whatever value of q; I also added the possibility of changing initial and final temperatures, and ah. However, if left off, the function will assume 500K and 350K respectively, as you had in your code:

data = Table[
         {q, function[q]},
         {q, {-0.1/60, -1/60, -10/60, -80/60.}}
       ];

ListPlot[
  data[[All, 2, 1]],
  PlotRange -> {{350, 500}, {400, 500}},
  ImageSize -> Medium, AxesLabel -> {"T", "Tf"}
]

ListPlot[
  data[[All, 2, 2]],
  PlotRange -> {{350, 500}, {-0.5, 3}},
  ImageSize -> Medium, AxesLabel -> {"T", "cp"}
]

plot of all Tf traces

plot of all cp traces

$\endgroup$
7
  • $\begingroup$ MarcoB the code is absolutely fantastic!! Thank you so much!. If I can ask one more question: How can I plot all the q's in the same plot (as to have 4 different graphs in the same plot)?. I am accepting your answer already and again I appreciate it very much, ! $\endgroup$
    – John
    Dec 6 '20 at 2:32
  • 1
    $\begingroup$ @John Sure, easiest will be to have the function return the values rather than the plots. See edit. $\endgroup$
    – MarcoB
    Dec 6 '20 at 2:44
  • $\begingroup$ MarcoB thank you so much!. I will have both versions of the code to consider. Again, thank you for your kind help ! $\endgroup$
    – John
    Dec 6 '20 at 2:47
  • $\begingroup$ @ MarcoB may I ask you another easy question for you. If I want to get the value of Tf[n] from the module to use it as the input for a second module. In this case, on heating I want to do exactly the same module that you constructed but instead of T0 being 350, I want T0 to be Tf[n] from the first module. How can I do this? $\endgroup$
    – John
    Dec 8 '20 at 0:40
  • 1
    $\begingroup$ @John Hmm maybe I'm not following. For each run, wouldn't Tf[n] be the last value of Tf? In other words, the function returns Tfc, which is a list of pairs of T[k] and Tf[k]. The last value of the last pair is Tf[n], right? So wouldn't function[yourQ][[1, -1, -1]] be Tf[n] for that run? $\endgroup$
    – MarcoB
    Dec 8 '20 at 1:46

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.