Message Boards Message Boards

How Do I combine NDSolve (RungeKutta4thOrder) and Loop(For or Do) together?

Posted 10 years ago

Hello everyone, previously I had ask about NDSolve with specified Runge Kutta 4th Order method for my equation. And here is what I got for combining few equations into one.

Z = 1

In[199]:= Clear[Z]

In[197]:= 
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]]


In[208]:= NDSolve[{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, {x, 0, 3}, 
 Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
   "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
 "StartingStepSize" -> 1/10]

Out[208]= {{y -> InterpolatingFunction[{{0., 3.}}, <>]}}

Z is a variable that I can change depending of the value for y[3] that I need. Here's the questions;

1) What command that I need to key-in in order for the equations to automatically generate the Plot and the value for y[3]? As what I have been currently doing is I keep on clicking "get solution", "plot" and " evaluate at point" buttons provided by Mathematica itself.
2) I need to use Loop command to find the value of Z so that I can get the value of y[3]=100 exactly. Currently I'm just try the value of Z randomly until I obtain y[3]=100 exactly. The For loop value of Z increment is 0.01 manually, the tutorials given to me was i++ command but the increment value is 1, how do I change it?

Thank you for your help guys.

POSTED BY: Thai Kee Gan
19 Replies

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.

POSTED BY: Ivan Morozov

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.

POSTED BY: Ivan Morozov

Also, you can modify f to return y instead of y[3.] (just remove [3.] in def. of f):

(* f[z][x] *)
f[1.5][3.]
f[#][3.] & /@ Range[1., 2., 0.1]
Plot[f[1.5][x], {x, 0, 3}]
ContourPlot[f[z][x], {x, 0, 3}, {z, 1, 3}]

I.M.

POSTED BY: Ivan Morozov

Hi,

you can add Twater as a parameter into def. of f :

(*Temperature of water*)
Twater = 300
(*Dimensions of the MHP*)
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 - 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
   ];
(* Tnew *)
170 + 2*z /. FindRoot[f[z,Twater] == 100., {z, 1}, Evaluated -> False]
(* 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]
(* Check *)
170 + 2*z /. FindRoot[f[z,Teq] == 100., {z, 1}, Evaluated -> False]

I.M.

POSTED BY: Ivan Morozov

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.

POSTED BY: Ivan Morozov

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.

POSTED BY: Ivan Morozov
Posted 10 years ago

Hello, I tried Break[] to get the value of z, when y[x] achieved 100,

For[z = 0,
  z < 100,
  z++,
  NDSolve[{Derivative[1][y][x] == y[x] + x^3 + 3 z, y[0] == 0}, 
   y, {x, 0, 3}]] If[y[3] > 100, Break[]]; y[3]

But the output is just y[3].

POSTED BY: Thai Kee Gan

To be sure that integration step is fixed you need to set "FixedStep" method as well as MaxStepFraction:

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,
    {x,0.,3.},
    z,
    Method -> {
       "FixedStep",
       Method -> {
         "ExplicitRungeKutta",
         "DifferenceOrder" -> 4, 
         "Coefficients" -> ClassicalRungeKuttaCoefficients
       }
    }, 
    StartingStepSize ->  #,
    MaxStepFraction -> #
] & ;

With this definition f is a function of step, z and x:

(* f[step][z][x] *)
f[0.5][1.][3.]
f[0.1][1.][3.]

You can then check that the step size is indeed fixed by requesting the grid in x used for integration:

f[0.5][1.]["Grid"]
f[0.1][1.]["Grid"]
(* commented output
{{0.}, {0.5}, {1.}, {1.}, {1.5}, {2.}, {2.}, {2.5}, {3.}, {3.}}    
{{0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6}, {0.7}, {0.8}, {0.9}, {1.}, {1.}, {1.1}, {1.2}, {1.3}, {1.4}, {1.5}, {1.6}, {1.7}, {1.8}, {1.9}, {2.}, {2.}, {2.1}, {2.2}, {2.3}, {2.4}, {2.5}, {2.6}, {2.7}, {2.8}, {2.9}, {3.}, {3.}}  *)

Also note that due to Piecewise[] there are almost identical points (auto discontinuity handling) on x-grid, like {1.}, {1.} (in fact {0.9999999999999716},{1.}).

I.M.

POSTED BY: Ivan Morozov
Posted 10 years ago

So, here's a little changes. For my equation, there's some information that I need to extract from Mathematica accordingly to the temperature of the water. For example, as the water temperature changes, the thermal conductivity of water will change as well .

I received error because Twater is not a value yet and I'm stuck.

Then again, sorry for the previous question that I did not explain fully for what I was trying to do.

 (*Temperature of water*)
     Twater = 300
     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] + Cwater/3 x^3 + 3 z ), 
         0 <= x <= 1}, {(y[x] + Cwater/3 x^2 + 2 z), 
         1 <= x <= 2}, {(y[x] + Cwater/3 x + z), 2 <= x <= 3}}], 
     y[0] == 0}, y[3.], {x, 0., 3.}, {z, Cwater}, 
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
      "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
    StartingStepSize -> 1/10];

 FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Tnew*)
 170 + 2*z /. FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*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]
 (*Check*)
 170 + 2*z /. FindRoot[f[z, Teq] == 100., {z, 1}, Evaluated -> False]
POSTED BY: Thai Kee Gan
Posted 10 years ago

Thank you, Ivan, for your help. This will definitely help me a lot in more than one way.

Is there a command that show the value of Z?, with the accuracy of my choice? for example I need to have the value of Z with the accuracy of 0.1 or 0.01 or 0.001, what command do I key-in?

By the way, I'm new to Mathematica, so I will never have think of this method.

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

POSTED BY: Thai Kee Gan
Posted 10 years ago

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.

POSTED BY: Thai Kee Gan
Posted 10 years ago

Hello Ivan,

I have encounter this error when I adjusted the starting point of my x0 value into higher or lower value. Is it not, when I have increase the starting value, I would not encounter the Maximum Steps error? I have tried few values and found the correct value that do not cause this error. But, even though I had this error message, the value of Z root can still be found.

FindRoot[f,{x,Subscript[x, 0]}] 

searches for a numerical root of f, starting from the point x=Subscript[x, 0].

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

ParametricNDSolveValue::mxst: Maximum number of 10000 steps reached at the point x$6275 == 0.6709750429824323`.

POSTED BY: Thai Kee Gan
Posted 10 years ago

Hello again,

If I have found the value of z

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]; FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

And Mathematica shows the answer in this form

{z -> 1.7049}

How do I use the value of the z to continue my equation? For example, the equation is A=100-z

POSTED BY: Thai Kee Gan
Posted 10 years ago

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]
POSTED BY: Thai Kee Gan
Posted 10 years ago

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?

POSTED BY: Thai Kee Gan
Posted 10 years ago

Thank you very much.

POSTED BY: Thai Kee Gan
Posted 10 years ago

Hello,

Please forgive me because I've tried some of the methods but it still doesn't work for me. Or, I just misunderstand some of the explanations

This is the original code that was written and I receive the error for Cwater = 0.61..... cannot be used as a parameter because Cwater gain a fixed value.

   (*Temperature of water*)
         Twater = 300
         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] + Cwater/3 x^3 + 3 z ), 
         0 <= x <= 1}, {(y[x] + Cwater/3 x^2 + 2 z), 
         1 <= x <= 2}, {(y[x] + Cwater/3 x + z), 2 <= x <= 3}}], 
     y[0] == 0}, y[3.], {x, 0., 3.}, {z, Cwater}, 
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
      "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
    StartingStepSize -> 1/10];

 FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Tnew*)
 170 + 2*z /. FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*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]
 (*Check*)
 170 + 2*z /. FindRoot[f[z, Teq] == 100., {z, 1}, Evaluated -> False]

Then I changed the command for Cwater by

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

And I received this error

ThermodynamicData::quants: List {t,300,360} does not consist of real numbers. >>

Shall I change into For loop with Temperature as i++ until 400 and stop when Tinitial = Tfinal? How do I do that?

POSTED BY: Thai Kee Gan
Posted 10 years ago

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.

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