Message Boards Message Boards

GROUPS:

Solve a system of differential equations in a loop

Posted 2 months ago
634 Views
|
7 Replies
|
0 Total Likes
|

Hello, What I try to do is. Solve a system of differential equations and obtain a set of values. Compare these values, in this case when it becomes >1, then change accordingly one initial value (with an If ) and repeat. Let's say I would like to do this five times. The code I use is the following. Αpparently, I do not implement the If command correctly in the code.

I hope this makes sense.

Thank you for your effort.

Best regards,

A.

PS. I have also attached the file.

 \[CapitalGamma] = 5/3;
 c4 = 6;
 smax = 50;
 Rinit = 1.05;
 \[Gamma]1 = 0.8;
 \[Gamma]2 = 0.9;


 t = Flatten@Table[i, {i, 1, 5, 1}];
 \[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
 Sol1 = Quiet[
     Catch[NDSolveValue[{D1[s] == 
         2 (-1 + R[s]) R[
           s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
            c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
             2 - \[CapitalGamma]) + 
            c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
              s]^-\[CapitalGamma]), 
        N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
            c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2)
              R[s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
               4 R[s] (-1 + \[CapitalGamma]) - 
               3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
             1 - \[CapitalGamma])), 
        R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
        Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
        R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
        WhenEvent[N1[s]* D1[s] == 0, Throw[{#, s, Abs[N1[s]/D1[s]] 
           , R[s]}]; "Stop Integration"]}, {R, \[Gamma]u}, {s, 0, 
        smax}]], NSDolveValue::noout] & /@ t
 j = Select[Sol1, #[[3]] > 1 &][[All, {1, 2, 3, 4}]];
 j[[All, 3]];
 If[j[[3, 3]] > 1, \[Gamma]2 = \[Gamma]uinit]
 If[j[[3, 3]] < 1, \[Gamma]1 = \[Gamma]uinit]
Attachments:
7 Replies

Agapi,

There is no need to loop. WhenEvent has the ability to change states when events happen and continue the integration. Look at the example in the documentation of WhenEvent. The first example is exactly what you need. I did not dive into your code to see what state you are checking in the If statement but I'll assume its [Gamma]u. You would do something like

WhenEvent[N1[s]* D1[s] == 0,If[\[Gamma]u > 1, \[Gamma]2 = \[Gamma]uinit,\[Gamma]1 = \[Gamma]uinit]]

Although I really do not understand what you are doing so you will have to change the action above to fit your situation. The documentation examples are good and you can execute multiple statements inside a list in the WhenEvent.

Regards

Neil

Posted 2 months ago

Thank you very much for your input. I have changed my code a bit using what you proposed and with the help of the documentations. What i can't do is to make the code start over a new integration using the new starting value of γuinit that meets the If condition.

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.81;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
{sol1, values} = 
  Reap@NDSolve[{D1[s] == 
      2 (-1 + R[s]) R[
        s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
          2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
          2 - \[CapitalGamma]) + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
          2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
           s]^-\[CapitalGamma]), 
     N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
         c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2)
           R[s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
            4 R[s] (-1 + \[CapitalGamma]) - 
            3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
          1 - \[CapitalGamma])), 
     R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
     WhenEvent[ {D1[s] == 0, N1[s] == 0} , 
      Sow[{s, D1[s], N1[s], Abs[N1[s]/D1[s]]}]; smax = s; 
      If[Abs[N1[s]/D1[s]] > 1, {\[Gamma]uinit = \[Gamma]2, 
        Print[\[Gamma]uinit]}]; 
      If[Abs[N1[s]/D1[s]] < 1, {\[Gamma]uinit = \[Gamma]1, 
        Print[\[Gamma]uinit]}]]}, {R, \[Gamma]u}, {s, 0, smax}];


   ParametricPlot[Evaluate[{R[s], \[Gamma]u[s]} /. sol1], {s, 0, smax}, 
 AspectRatio -> 1, PlotStyle -> ColorData["GreenPinkTones"], 
 AxesLabel -> {"R", "\[Gamma]u"}]

Your bug was that you included the line smax = s in the WhenEvent[]. By doing this you forced the integration to stop just when you wanted the If statement to alter your constants.

This works fine:

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.81;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
{sol1, values} = 
  Reap@NDSolve[{D1[s] == 
      2 (-1 + R[s]) R[
        s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(2 - \[CapitalGamma]) + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^-\[CapitalGamma]), 
     N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
         c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2) R[
            s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
            4 R[s] (-1 + \[CapitalGamma]) - 
            3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(1 - \[CapitalGamma])), 
     R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
     WhenEvent[{D1[s] == 0, 
       N1[s] == 0}, {Sow[{s, D1[s], N1[s], Abs[N1[s]/D1[s]]}],
       If[Abs[N1[s]/D1[s]] > 1, \[Gamma]uinit = \[Gamma]2; 
        Print[\[Gamma]uinit]];
       If[Abs[N1[s]/D1[s]] < 1, \[Gamma]uinit = \[Gamma]1; 
        Print[\[Gamma]uinit]]}]}, {R, \[Gamma]u}, {s, 0, smax}];

Regards,

Neil

Posted 2 months ago

Thank you again for your help, but this is not doing what I was expecting. I think I did not explain my goal correctly. I want the integration to stop the time the condition D[1]==0,N1[s]==0 is met. This value i name smax. Then it should go back at the beginning with a new initial value of γuinit and start the integration all over (s=0) again. The visual result should be a different line close to the original that it will stop at a new smax, when the same condition is met. Rinse and repeat let's say n times.

Its going to be something like this: (if I understand what you want)

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.81;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
sols = Table[
  Reap@NDSolve[
    smax = 2; {D1[s] == 
      2 (-1 + R[s]) R[
        s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(2 - \[CapitalGamma]) + 
         c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
             2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^-\[CapitalGamma]), 
     N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
         c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2) R[
            s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
            4 R[s] (-1 + \[CapitalGamma]) - 
            3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^(1 - \[CapitalGamma])), 
     R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
     R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit, 
     WhenEvent[{D1[s] == 0, N1[s] == 0}, 
      Sow[{s, D1[s], N1[s], Abs[N1[s]/D1[s]]}]; smax = s;
      If[Abs[N1[s]/D1[s]] > 1, {\[Gamma]uinit = \[Gamma]2, 
        Print[\[Gamma]uinit]}];
      If[Abs[N1[s]/D1[s]] < 1, {\[Gamma]uinit = \[Gamma]1, 
        Print[\[Gamma]uinit]}]]}, {R, \[Gamma]u}, {s, 0, smax}], 
  5]; 
data =  Map[({R[s], \[Gamma]u[s]} /. #) &, 
      sols[[All, 1, 1]]]; 
ParametricPlot[data, {s, 0, smax}, 
     AspectRatio -> 1, AxesLabel -> {"R", "\[Gamma]u"}]

The idea is to make a list of solutions and then plot the list of results.

Regards,

Neil

Posted 2 months ago

Τhank you for the help Neil.

The list and the plot works as it should. The "If" doesn't do what I want it to do. The idea is to change the γuinit to (γ1+(γ1+γ2)/2)/2 inside the first If or to (γ2+(γ1+γ2)/2)/2 inside the second, in order to bisect the initial interval of γ1,γ2. I tried the Piecewise also instead of if ,after reading some posts, but did not work either. I used the last line in order to visual see the results, but this only shows the result of the last integration. Is it possible to see it for the previous integrations? Any ideas would be really appreciated. Regards, Agapi

\[CapitalGamma] = 5/3;
c4 = 6;
smax = 2;
Rinit = 1.05;
\[Gamma]1 = 0.8;
\[Gamma]2 = 0.9;
\[Gamma]uinit = (\[Gamma]1 + \[Gamma]2)/2;
sol1 = Table[
   Reap@NDSolve[{D1[s] == 
       2 (-1 + R[s]) R[
         s] (-(-1 + \[CapitalGamma]) \[Gamma]u[s] + 
          c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
           2) (-2 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
           2 - \[CapitalGamma]) + 
          c4 ((-1 + R[s]) R[s]^3)^((1 - \[CapitalGamma])/
           2) (-1 + \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[
            s]^-\[CapitalGamma]), 
      N1[s] == -(1 + \[Gamma]u[s]^2) (1 - \[CapitalGamma] + 
          c4 (-1 + R[s])^(1/2 - \[CapitalGamma]/2)
            R[s]^(-(3/2) (-1 + \[CapitalGamma])) (2 + 
             4 R[s] (-1 + \[CapitalGamma]) - 
             3 \[CapitalGamma]) \[CapitalGamma] \[Gamma]u[s]^(
           1 - \[CapitalGamma])), 
      R'[s] == D1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
      Derivative[1][\[Gamma]u][s] == N1[s]/Sqrt[D1[s]^2 + N1[s]^2], 
      R[0] == Rinit, \[Gamma]u[0] == \[Gamma]uinit,
      WhenEvent[ {D1[s] == 0, N1[s] == 0} , 
       j = Sow[{s, N1[s], D1[s], Abs[N1[s]/D1[s]]}]; smax = s;
       If[
        Abs[N1[s]/D1[s]] > 
         1, {\[Gamma]uinit = (\[Gamma]1 + \[Gamma]uinit)/2, 
         Print[\[Gamma]uinit]}];
       If[
        Abs[N1[s]/D1[s]] < 
         1, {\[Gamma]uinit = (\[Gamma]2 + \[Gamma]uinit)/2, 
         Print[\[Gamma]uinit]}]]}, {R, \[Gamma]u}, {s, 0, smax}], 5];
data = Map[({R[s], \[Gamma]u[s]} /. #) &,
   sol1[[All, 1, 1]]];
ParametricPlot[data, {s, 0, smax}, AspectRatio -> 1, 
 AxesLabel -> {"R", "\[Gamma]u"}]
j[[All]] & /@ sol1

I really don't understand what you are trying to do. Can you give a simple example? At first I thought you wanted to hit an event and then change some settings and continue. Next you said that you only wanted to do one simulation and stop at the event, reset everything, change a setting and start over. This is what your code is currently doing. Four of the simulations are identical because once you make your setting change the answer is always the same. Your code runs in both cases and produces results so without really understanding what you want, it's hard to be more helpful.

Regards

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