# Randomizing Initial Conditions of NDSolve Differential Equation Phase Space Plots with Looping Structures

I am looking for a way to randomize the initial conditions on the numerical approximation for the angular displacement of a plane pendulum. I have been trying for hours, and believe a For loop is easiest, but I don't know how to get it to work.

My objective is to produce 10 total phase plots, overlaid on each other, of different initial conditions with 5 of the plots of θ_0 between (-π and π) and 5 plots of θ_0 outside of (-π and π).

When I run this, I get a list of plot[11] and then some broken/wierd stuff that aren't nicely overlaid on plot[11].

### Edit

It appears an ic2a/b value of +/- 5 is too large. Sometimes the initial conditions break the orbits.

### Edit 2

I forgot to define g and R. Fixed.

### Edit 3

Fixed the code. See comments below for diagnosis and steps to solution. Thanks for the help. RandomColor still doesn't manifest in the plots, but it's superficial (as is a legend) in this situation. Also I haven't had a phase plot to be shifted from the origin along the [Theta] axis. I still consider the issue solved because I have my plots and am at least able to generate random values for the initial conditions.

g=9.8;
R=1;
SeedRandom[RandomInteger[100]]; (*sets a random seed to generate below RandomReal's*)

Array[plot, 11];

plot[11] =
Plot[f[x], {x, -100, 100},
PlotRange -> {{-20, 20}, {-20, 20}},
Frame -> True,
FrameLabel ->
{"Displacement, θ", Row[{"Momentum, ", Subscript[P, θ]}]},
PlotLabel ->
"Phase Space of Approoximate Motion of A Vertical Plane Pendelum",
ImageSize -> Large];

For[i = 1, i < 11, i++,
If[i < 6,
ic1a = RandomReal[{-π, π}];
(*5 initial condition displacements between-π and π*)
ic2a = RandomReal[{-5, 5}];
(*initial condition momentum=random number between-5 and 5*)
theta =
NDSolve[
{x''[t] == -(g/R) Sin[x[t]], x[0] == ic1a, x'[0] == ic2a},
x, {t, 0, 15},
MaxSteps -> 20000];
plot[i] =
ParametricPlot[Evaluate[{x[t], Derivative[1][x][t]} /. theta], {t, 0, 10},
AspectRatio -> 1,
PlotLabel -> Style["Phase Space"],
AxesLabel -> {"θ", "θ'"},
PlotStyle -> RandomColor];,
(*else,i.e.i≥6*)
ic1b = RandomReal[{π, -π}];
(*5 initial condition displacements outside of-π and π*)
ic2b = RandomReal[{-5, 5}];
(*initial condition momentum=random number between-5 and 5*)
theta =
NDSolve[
{x''[t] == -(g/R) Sin[x[t]], x[0] == -(ic1b), x'[0] == ic2b},
x, {t, 0, 15},
MaxSteps -> 20000];
plot[i] =
ParametricPlot[Evaluate[{x[t], Derivative[1][x][t]} /. theta], {t, 0, 10},
AspectRatio -> 1,
PlotLabel -> Style["Phase Space"],
AxesLabel -> {"θ", "θ'"},
PlotStyle -> RandomColor];]]

Show[plot[11], Table[plot[i], {i, 1, 10}]]

• Did your edit solve your problem? BTW, I tried running your code and ran into many undefined parameters / plots. If you still have a question, please edit your code so that people can run it. – Chris K Oct 17 '18 at 11:55
• ps. Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! – Chris K Oct 17 '18 at 11:55
• Hi Josh, your use of BlockRandom around a RandomReal will basically reseed the generator, which means you will get the same result for all samples if you call it before every use of RandomReal, which is probably not what you want. Try to call SeedRandom with a fixed seed once at the start of your program (outside the For loop) for reproduceability and then inside your loop call RandomReal freely to get different random samples. – Thies Heidecke Oct 17 '18 at 14:57
• Oops, I forgot to define g and R @Chrisk. Let g = 9.8 and R =1. I've updated in the code. – Josh Albert Oct 18 '18 at 0:45
• @ThiesHeidecke that worked! Thanks so much. That and defining fixed the generation of plots. Now I'm just having trouble showing all the plots on the same base plot, plot[11]. – Josh Albert Oct 18 '18 at 0:51

Because i think this will make your life easier with Mathematica in the future i did a few refactorings on your solution to make it more idiomatic and give you some ideas what simpler options you have sometimes to get your result. First the fully refactored code:

ODESolFromICs[{x0_, xd0_}] := NDSolve[
{x''[t] == -(g/R) Sin[x[t]], x[0] == x0, x'[0] == xd0}
, x, {t, 0, 15}, MaxSteps -> 20000
];

PlotSolution[sol_, col_] := ParametricPlot[
Evaluate[{x[t], Derivative[1][x][t]} /. sol]
, {t, 0, 10}, PlotStyle -> col
];

plotoptions = {
PlotRange -> {{-20, 20}, {-10, 10}},
Frame -> True,
FrameLabel -> {"Displacement, \[Theta]", Row[{"Momentum, ", Subscript[P, \[Theta]]}]},
PlotLabel -> "Phase Space of Approoximate Motion of A Vertical Plane Pendelum",
ImageSize -> Large,
AspectRatio -> Automatic
};

g = 9.8;
R = 1;
SeedRandom[1];(*choose a different number for a different randomized run*)

plots = Table[
initialangle = If[i < 6,
RandomReal[{-\[Pi], \[Pi]}](*5 initial condition displacements between-\[Pi] and \[Pi]*),
First@RandomPoint[MeshRegion[{{-5}, {-\[Pi]}, {\[Pi]}, {5}}, Line[{1, 2}],Line[{3, 4}]}]](*5 initial condition displacements outside of-\[Pi] and \[Pi]*)
];
initialmomentum = RandomReal[{-5, 5}];(*initial condition momentum=random number between-5 and 5*)
theta = ODESolFromICs[{initialangle, initialmomentum}];
PlotSolution[theta, RandomColor[]]
, {i, 10}
];

Show[plots, plotoptions]


Now for what is changed:

• First we factor out the ODE solving and plotting the solution into their own functions ODESolFromICs and PlotSolution. This reduces code duplication and makes the code more readable. Also we added a parameter to give a plot color to ParametricPlot. Same holds true for the initial conditions for NDSolve.
• Array[plot,11]; is in general not needed as 'initialization', as the meaning of the expression plot[5] will be given by the plot[5]=... (see Set in the documentation) replacement rule. We don't have to declare those expression before we use them, essentially because they are recognized patterns instead of uninitialized variables. Also the semicolon will discard the result, which means the line doesn't have any lasting effect.
• We can replace the For loop by a Table which iterates i from 1 to 10 and save the resulting list of plots. This replaces the previous different plot[i] expressions.
• We don't need a plot[11] expression that just holds the plot options we want, instead we can save those in a list and apply them at the end when we call Show. Directly giving those options to Show would also work.
• We can move in the If[i<6,...] to just influence the initial conditions, which makes the code easier to understand

Hope this helps! :)