Hello again,

Now I have tried this method manually, it's suppose to work. Just wondering how do I do it automatically. For example

(*Temperature of water*)
Twater = 500
Cwater = QuantityMagnitude[
  ThermodynamicData["Water", 
   "ThermalConductivity", {"Temperature" -> 
     Quantity[Twater, "Kelvins"]}]]

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 - 18000*Cwater), 
        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.};
Show[ListPlot[Evaluate[{#, f[#]} & /@ Range[60., 100., 0.05]], 
  Joined -> True, Frame -> True, FrameLabel -> {"Z", "Y[3]"}], 
 Graphics[{PointSize[Large], Red, Point[point]}], ImageSize -> 400]

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

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

How do I do this command?:

Tinitial. Found Tnew. Tinitial = Tnew, False. Replace Tinitial with Tnew Repeat. Tinitial = Tnew, True. Stop. OR Stop when iteration is 100 and no suitable value is found

NOTE: Tried this equation with few iterations, Tnew keep getting further away from Tinitial.

You can define Cwater as a function of temp. :

    Cwater = QuantityMagnitude[
        ThermodynamicData[
           "Water",
           "ThermalConductivity", {"Temperature" ->Quantity[#, "Kelvins"]}
        ]
    ] & ;
    (* eg: *)
    Cwater[300.]
    Cwater[310.]
    Cwater[320.]

Then f : {z,c} -> y[3] where c is a parameter (former Cwater) can be called as :

f[1.9,Cwater[300.]]

I.M.

Btw, what does # represents?

# is a slot, it is used to define pure functions, e.g. :

f1[x_] = x^2 ;
f2 = #^2 & ;
(* the returned result is the same *)
{f1[2],f2[2]}
(* internal representation is different *)
DownValues[f1]
DownValues[f2]

See this post for more info as well as documentation.

Here is a cleaner (I hope) version of your code:

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]
] ;

(* Define function f : {z,t} ->  y[3.] *)
f = ParametricNDSolveValue[
    {
    Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - t), 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,t},
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];

(* Define function g : t_initial -> t_final *)
g = 170 + 2 z /. FindRoot[f[z,#] == 100., {z, 1}, Evaluated -> False] & ;

(* Find fix point of g, i.e. g[T] = T *)
Tinitial = 300. ;
Tfinal = x /. FindRoot[g[x] == x,{x,Tinitial},Evaluated -> False]
(* Or *)
Quiet[FixedPoint[g,Tinitial]]

I.M.

Thanks again, Ivan.

May I know how does the equation work? At first, the equation help me solve z with initial Twater = 300

f = ParametricNDSolveValue[{Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - t), 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,t},
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];

Then it help me to look for a Tnew

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

With the value is not True, The third function repeat ParametricNDSolveValue until it reaches True

(* Find eq. temp. *)
g = 170 + 2 z /. FindRoot[f[z,#] == 100., {z, 1}, Evaluated -> False] & ;
Teq = x /. FindRoot[g[x]==x,{x,Twater},Evaluated -> False]

Then finally, as it reaches True, Recheck the Twater and conclude the correct value for Twater= Tnew is True

(* Check *) 170 + 2*z /. FindRoot[f[z,Teq] == 100., {z, 1}, Evaluated -> False]

Is that how it works? Btw, what does # represents?

Hello again,

For this equation, I would like to do an iteration to increase Twater temperature until the Tnew are the same. I would like to know which type of loop that I can use until Twater=Tnew Shall I use ParametricNDSolve or other type of method?

(*Temperature of water*)
Twater = 300

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 - Twater), 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.};
Show[ListPlot[Evaluate[{#, f[#]} & /@ Range[60., 100., 0.05]], 
  Joined -> True, Frame -> True, FrameLabel -> {"Z", "Y[3]"}], 
 Graphics[{PointSize[Large], Red, Point[point]}], ImageSize -> 400]

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

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

Thanks again Ivan. Your guidance really help me a lot. I'm a student taking masters right now, so I need to confirm with the methods that I am using to solve the equations. For my method of solution, I have determined to use RK4 with certain fixed step in my work, to achieve high accuracy.

Just to make confirmation, is ParametricNDSolveValue[] command, also use Runge-Kutta 4th Order Method as the base of the solution?

No, but setting method is identical to that of NDSolve[], also do you need a fixed step rk4?

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),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
];

Is there a command that show the value of Z?

The defined f functions takes value of Z and returns y[3] as output. Than the code:

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

finds the value of Z for which y[3] = 100.

This creates a list {{z1,y3[z1]},{z2,y3[z2]},...} (see documentation for Range[] function)

{#, f[#]} & /@ Range[1., 3., 0.05]

Note, f can be plotted as a regular function of Z:

Plot[f[z], {z, 1, 3}]

I.M.

Hi Thai,

Instead of loops try ParametricNDSolve[].

f = ParametricNDSolveValue[
    {
       Derivative[1][y][x] == 
       Piecewise[
         {
          {(y[x] + x^3 + 3 z),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
] ;

point = {z /. FindRoot[f[z] == 100.,{z,1},Evaluated->False], 100.} ;
Show[
    ListPlot[Evaluate[{#,f[#]} &/@ Range[1.,3.,0.05]],Joined->True,Frame->True,FrameLabel->{"Z","Y[3]"}],
    Graphics[{PointSize[Large],Red,Point[point]}],
    ImageSize -> 400
]

I.M.