# Two different variables in a Do Loop

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?

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"}
]


• 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, !
– John
Dec 6 '20 at 2:32
• @John Sure, easiest will be to have the function return the values rather than the plots. See edit. Dec 6 '20 at 2:44
• MarcoB thank you so much!. I will have both versions of the code to consider. Again, thank you for your kind help !
– John
Dec 6 '20 at 2:47
• @ 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?
– John
Dec 8 '20 at 0:40
• @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? Dec 8 '20 at 1:46