# Solve a system of differential equations in a loop

Posted 8 months ago
1365 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
Sort By:
Posted 8 months ago
 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.RegardsNeil
Posted 8 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"}] 
Posted 8 months ago
 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 8 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.
Posted 8 months ago
 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
 Τ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