Message Boards Message Boards

How do I repeat/loop multiple functions that are related?

Posted 9 years ago

Hello guys, I'm stuck with this question for some time. For my equation, I need to use WolframAlpha for my data. Firstly, I fix my initial water temperature as 500 Kelvin. Eg:

Twater = 300
Cwater = QuantityMagnitude[
  ThermodynamicData["Water", 
   "ThermalConductivity", {"Temperature" -> 
     Quantity[Twater, "Kelvins"]}]]

With Cwater determined, I can implement the next 2 functions

Dwater = (5 Cwater)/\[Pi]
Ewater = Dwater^3 + 2 Cwater

Then by using ParametricNDSolveValue, I can determine the value of Z

ClassicalRungeKuttaCoefficients[4, prec_] := 
  With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
    bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
   N[{amat, bvec, cvec}, prec]];

f = ParametricNDSolveValue[{Derivative[1][y][x] == 
     Piecewise[{{(y[x] + x^3 + 3 z - 120*Ewater), 0 <= x <= 1},
       {(y[x] + x^2 + 2 z), 1 <= x <= 2},
       {(y[x] + x + z), 2 <= x <= 3}}],
    y[0] == 0},
   y[3.],
   {x, 0., 3.},
   z,
   Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
   StartingStepSize -> 1/10];

point = {z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False], 
   100.};

FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

Lastly, to find the new temperature of water

Tnew = 170 + 1.89*z /. 
  FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

I wanted to repeat Twater with Tnew until Twater=Tnew->True and the loop stop. I need to start Twater with a value, (only to replace Twater=Tnew after the first equation and keep press shift+enter repeatedly)

I've tried some Do Loop or For Loop and even FixedPoint but still don't know how to combine multiple functions/equations into 1 loop.... I apologize because I've been asking this question a few times, or maybe some of you have answered me in some way, but I still did not understand how to do this. Thank you very much for your time.

POSTED BY: Thai Kee Gan
12 Replies

Mathematica has a short tutorial entitled "Arbitrary-Precision Numbers." You can find it by searching for "tutorial/ArbitraryPrecisionNumbers" in Help/Wolfram Documentation or by typing SetPrecision, highlighting it, pressing F1, then find the link for the tutorial at the bottom of the page.

POSTED BY: Charles Elliott
Posted 9 years ago

Thanks again for your help and your advice, Dr. Charles Elliott (correct me if I'm wrong). I will try to digest everything. Hopefully.

POSTED BY: Thai Kee Gan

I made prog run, but it still does not converge. I believe that the coefficients to the ExplicitRungeKutta method may be incorrect. The four changes I made are:

  1. There was no initial value for point, so the inner While never executed.

  2. Since the RHS of ClassicalRungeKuttaCoefficients[4, prec] was never executed, the coefficients had no value. ClassicalRungeKuttaCoefficients[4, prec] is a non-starter; you cannot define a constant in the definition of a function like that. In any case, it was not being called as a function in NDSolve, so the coefficients were not defined. "point = y[1] /. f[[1]];" was doing exactly what it was supposed to; however, f had no meaningful value;

  3. The lack of a ';' after the definition of NDSolve, meant that f was interpreted as a multiplication of the Print stmt after it; so, f had the value garbage * null, null being what Print returns.

  4. As you sent it, prog would not define because of an imbalance in [ ]. I made several changes to balance the brackets. You might want to check that the logic still makes sense.

Recommendations:

You may be better off using the AccuracyGoal and/or PrecisionGoal options to NDSolve before using SetPrecision there. In any case, SetPrecison conflicts with the "MinPrecison = 50;" stmt in the beginning of the definition of prog. In theory, $MinPrecison controls the number of digits of all computations in the Module. Mathematica has a short tutorial (entitled, "Arbitrary-Precision Numbers") on the use of SetPrecision and SetAccuracy built-in. Just highlight one or the other and press F1. A link to the tutorial is at the bottom.

I put in code to save and restore the value of $MinPrecision.

When you add Print stmts to the routine, put in a label, such as Print["Point: ", point, "] instead of Print[point].

When you define a While or other kind of loop, try to remember to put the closing "];" on its own line and perhaps even comment it, thusly:

]; (* End inner while *)

That way, Mathematica will indent the lines properly for you so you can clearly see where the loops start and end. It is almost impossible to put in too many comments. A month from now, you will not remember what you did if you don't put in comments. Consider you are building a library of working code.

I suggest you make NDSolve and the RungeKutta method work outside of the prog Module. The edit prog, and try again.

Instead of defining a function, such as ClassicalRungeKuttaCoefficients[prec_Integer] := ..., within a module consider using anonymous functions instead:

Function[prec, body]

Then put the Function[prec, body] within the NDSolve. Anonymous functions are much, much faster in execution. But that is frosting on the cake. Make it work first, then make it faster.

Attachments:
POSTED BY: Charles Elliott
Posted 9 years ago

Hello again, sorry for the trouble. I've tried NDSolve method to search for "z" instead of ParametricNDSolveValue So here's what I done but still there's some error.

prog := Module[{Twater = 300., Tnew = 0, Cwater, Dwater, Ewater, 
    ClassicalRungeKuttaCoefficients, f, x, y, z = 60, point, ans, 
    tol = 0.001, iters = 0, iterLim = 100, prevTnew = 10^5, 
    plotList = {}}, $MinPrecison = 50;
   Cwater = 
    QuantityMagnitude[
     ThermodynamicData["Water", 
      "ThermalConductivity", {"Temperature" -> 
        Quantity[Twater, "Kelvins"]}]];
   ans = "Twinkle, twinkle, little star,
             How I wonder what you are.
             Up above the world so high,
             Like a diamond in the sky.
             Twinkle, twinkle, little star,
             How I wonder what you are!";

   While[(Abs[prevTnew - Tnew] > tol) && (iters++ < iterLim), 
    Dwater = (5 Cwater)/Pi;
    Ewater = Dwater^3 + 2*Cwater;

    While[
      (Abs[point] < 100) && (iters++ < iterLim),

      ClassicalRungeKuttaCoefficients[4, prec_] :=
       With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
         bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
        N[{amat, bvec, cvec}, prec]];

      f = 
       NDSolve[{SetPrecision[
          Derivative[1][y][x] == 
           Piecewise[{{(y[x] + x^3 + 3 z - 120*Ewater), 
              0 <= x <= 1}, {(y[x] + x^2 + 2 z), 
              1 <= x <= 2}, {(y[x] + x + z), 2 <= x <= 3}}], 40], 
         y[0] == 0}, y, {x, 0., 3.}, 
        Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
          "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
        StartingStepSize -> 1/100];

      point = y[1] /. f[[1]];

      z++

      ]

     Print[point];



    prevTnew = Tnew;

    Tnew = 170 + 1.89*z;

    Cwater = 
     QuantityMagnitude[
      ThermodynamicData["Water", 
       "ThermalConductivity", {"Temperature" -> 
         Quantity[Tnew, "Kelvins"]}]];

    (*Note Tnew here!*)
    Print["Twater: ", Twater, " Tnew: ", Tnew];];

   (*End while*)

   ans = {"Tnew: ", Tnew, "Cwater: ", Cwater, "Dwater: ", Dwater, 
     "Ewater: ", Ewater};
   Print["The answer is:\n", ans];];

What I tried to do is to have "While Loop" inside a "While Loop" and instead of FindRoot, I use ReplaceAll, Here's the plan,

The inner loop is to find the value of "z" so that the "point" is 100, and the increment of "z" is 1 each time. After it reach the correct value of "z", the the command will continue to look for Tnew, and repeat again.

But the problem is

 point = y[1] /. f[[1]];

does not work properly in the While loop. I can do it outside of While loop but not inside. So, i'm stuck here

POSTED BY: Thai Kee Gan
POSTED BY: Charles Elliott
Posted 9 years ago

Hello Charles, thank you for your advice. You are correct about my background, that i have not been exposed to numerical analysis as much as I supposed to have. In the future, I will definitely take up PhD program if there is a chance.

POSTED BY: Thai Kee Gan
Posted 9 years ago

OMG THANK YOU VERY MUCH. It totally works. I will not be able to get it correct by myself. The codes posted here is an alternate code for my research. Sorry I cant disclose it here. Here's the explanation. I'm an engineering master student working on heat transfer. In this case, it's about fluid. For the code, I would need to find the value of "z" so that "y=100". For a given condition, the "Twater" changes and the thermal conductivity of water. Then, "z" will be replace into Tnew equation. If Twater=Tnew, then that is the temperature that we want. At first, I was thinking to try to iterate with If,Do,While loop, with step size 0.001...., but that will take forever. After some trial, Tnew will actually converge, so I wanted to replace Twater with Tnew. This equation is some kind of temperature loop.

I will study the code given and learn more about the commands.

BTW,

Cwater = QuantityMagnitude[
   ThermodynamicData["Water", 
    "ThermalConductivity", {"Temperature" -> 
      Quantity[Twater, "Kelvins"]}]];

Is there any reason to have Cwater outside of the While loop? Well I tried putting it into the While loop and the loop just stop after 2 iterations. It will be helpful for me to understand this.

POSTED BY: Thai Kee Gan

I made the algorithm converge to an answer; I have no idea if it is correct. Here are the changes I made:

I replaced the statement Twater = Tnew, with Cwater = QuantityMagnitude[ ThermodynamicData["Water", "ThermalConductivity", {"Temperature" -> Quantity[Tnew, "Kelvins"]}]]; <-- note the Tnew here. I added the variable prevTnew, initialized it to 10^5, and set it to Tnew just before Tnew is computed. I changed the completion test to (Abs[prevTnew - Tnew] > tol) && (iters++ < iterLim). Also, please note I created the variable plotList, and initialized it to plotList = {}. Then after you compute point, I added AppendTo[plotList, point]. Later, you could put in Print[ListPlot[plotList]] if the reason you compute point is to plot it. I don't understand the 100 in point. Here is the new routine:

prog := Module[{Twater = 300., Tnew = 0, Cwater, Dwater, Ewater, 
    ClassicalRungeKuttaCoefficients, f, x, y, point, ans, tol = 0.001,
     iters = 0, iterLim = 100, prevTnew = 10^5, plotList = {}},
   $MinPrecison = 50; 
   Cwater = 
    QuantityMagnitude[
     ThermodynamicData["Water", 
      "ThermalConductivity", {"Temperature" -> 
        Quantity[Twater, "Kelvins"]}]];
   ans = "Twinkle, twinkle, little star,
         How I wonder what you are.
         Up above the world so high,
         Like a diamond in the sky.
         Twinkle, twinkle, little star,
         How I wonder what you are!";
   While[(Abs[prevTnew - Tnew] > tol) && (iters++ < iterLim), 
    Dwater = (5 Cwater)/Pi;
    Ewater = Dwater^3 + 2*Cwater;
    ClassicalRungeKuttaCoefficients[4, prec_] := 
     With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
       bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
      N[{amat, bvec, cvec}, prec]];
    f = ParametricNDSolveValue[{SetPrecision[
        Derivative[1][y][x] == 
         Piecewise[{{(y[x] + x^3 + 3 z - 120*Ewater), 
            0 <= x <= 1}, {(y[x] + x^2 + 2 z), 
            1 <= x <= 2}, {(y[x] + x + z), 2 <= x <= 3}}], 40], 
       y[0] == 0}, y[3.], {x, 0., 3.}, z, 
      Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
        "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
      StartingStepSize -> 1/10];
    point = {z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False], 
      100.};
    AppendTo[plotList, point];
    Print[point];
    FindRoot[f[z] == 100., {z, 1}, Evaluated -> False];
    prevTnew = Tnew;
    Tnew = 
     170 + 1.89*z /. 
      FindRoot[f[z] == 100., {z, 1}, Evaluated -> False];
    Cwater = 
     QuantityMagnitude[
      ThermodynamicData["Water", 
       "ThermalConductivity", {"Temperature" -> 
         Quantity[Tnew, "Kelvins"]}]]; (* Note Tnew here! *)
    Print["Twater: ", Twater, " Tnew: ", Tnew];
    ];(*End while*)
   ans = {"Tnew: ", Tnew, "Cwater: ", Cwater, "Dwater: ", Dwater, 
     "Ewater: ", Ewater};
   Print["The answer is:\n", ans];
   ];
POSTED BY: Charles Elliott
Posted 9 years ago

Hold on..... I suppose when we replace Twater with Tnew, then

While[(Abs[Twater - Tnew] > tol) && (iters++ < iterLim),

will not hold True anymore, that is why the iteration only run once.... So, how do I change it?, or I can just use the conventional method, which is much more "messy"?

n = 5; While[0 < n - 1 < 9, Print[n]; n = n - 0.5]
POSTED BY: Thai Kee Gan
Posted 9 years ago
POSTED BY: Thai Kee Gan

The code below is the general method, but I am not sure it achieves the correct answer. Start it by "Shift+Enter" or "Keypad Enter" on the cell containing prog. You debug it by putting in print statements until you find the error.

prog := Module[{Twater = 300, Tnew = 0, Cwater, Dwater, Ewater, 
    ClassicalRungeKuttaCoefficients, f, x, y, point, ans, tol = 0.001,
     iters = 0, iterLim = 100},
   Cwater = 
    QuantityMagnitude[
     ThermodynamicData["Water", 
      "ThermalConductivity", {"Temperature" -> 
        Quantity[Twater, "Kelvins"]}]];
   ans = "Twinkle, twinkle, little star,
     How I wonder what you are.
     Up above the world so high,
     Like a diamond in the sky.
     Twinkle, twinkle, little star,
     How I wonder what you are!";
   While[(Abs[Twater - Tnew] > tol) && (iters++ < iterLim),
    Dwater = (5 Cwater)/Pi;
    Ewater = Dwater^3 + 2*Cwater;
    ClassicalRungeKuttaCoefficients[4, prec_] := 
     With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
       bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
      N[{amat, bvec, cvec}, prec]];
    f = ParametricNDSolveValue[{Derivative[1][y][x] == 
        Piecewise[{{(y[x] + x^3 + 3 z - 120*Ewater), 
           0 <= x <= 1}, {(y[x] + x^2 + 2 z), 
           1 <= x <= 2}, {(y[x] + x + z), 2 <= x <= 3}}], y[0] == 0}, 
      y[3.], {x, 0., 3.}, z, 
      Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
        "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
      StartingStepSize -> 1/10];
    point = {z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False], 
      100.};
    Print[point];
    FindRoot[f[z] == 100., {z, 1}, Evaluated -> False];
    Tnew = 
     170 + 1.89*z /. 
      FindRoot[f[z] == 100., {z, 1}, Evaluated -> False];
    Print["Twater: ", Twater, " Tnew: ", Tnew];
    Twater = Tnew;
    ];  (* End while *)
   Print["The answer is:\n", ans];
   ];

prog

POSTED BY: Charles Elliott
Posted 9 years ago

I was wondering if Module[] has something to do with it. It's just that the examples given in the documentation was just for simple tasks. I'm not in my Uni's lab right now (since the software was bought by Uni), but I'll get to it as soon as I can and let you know the result. Thank you very much.

POSTED BY: Thai Kee Gan
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract