Group Abstract Group Abstract

Message Boards Message Boards

Solve a system of differential equations in a loop

Agapi Pop
Posted 5 years ago

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:
NLoop3.nb
POSTED BY: Agapi Pop
7 Replies
Sort By:

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
POSTED BY: Neil Singer

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 BY: Neil Singer
Agapi Pop
Agapi Pop
Posted 5 years ago
POSTED BY: Agapi Pop

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 BY: Neil Singer
Agapi Pop
Agapi Pop
Posted 5 years 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.
POSTED BY: Agapi Pop

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 BY: Neil Singer
Agapi Pop
Agapi Pop
Posted 5 years 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"}]
POSTED BY: Agapi Pop
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard