# Running a Numerical Simulation on a System of Differential Equations with unique initial conditions

I have a system of differential equations as follows:

β*M*L^2*D[Θ[i,t],{t,2}]=M*g*μ*L*(-Sin[Θ[i,t]])+τ[Θ[i+1,t],Θ[i,t]], (*for i=1*)
β*M*L^2*D[Θ[i,t],{t,2}]=M*g*μ*L*(-Sin[Θ[i,t]])+τ[Θ[i-1,t],Θ[i,t]]+τ[Θ[i+1,t],Θ[i,t]], (*for i between 1 and n-1*)
β*M*L^2*D[Θ[i,t],{t,2}]=M*g*μ*L*(-Sin[Θ[i,t]])+τ[Θ[i-1,t],Θ[i,t]], (*for i=n*)


Basically, what the equation describes is a series of coupled oscillators from 0 to n, with the i^th oscillator being connected to the (i-1)^th and (i+1)^th ones with a spring. τ[Θ[j],Θ[i]] describes the torque exerted by the j^th oscillator on the i^th.

Edit: My initial conditions for the systems are as follows.

• Initially, at time t=0, all the oscillators have:
• an angular position of 0, that is, Θ[i,0]=0 for all i from 0 to n
• and an angular velocity of 0, that is D[Θ[i,t],{t,1}]=0 at t=0 for all i from 0 to n
• Now, what I do next is to apply a force on the first oscillator causing it to accelerate and follow a certain path of motion over a specific initial time interval. For example, the force that I apply could cause it to complete to full rounds from time t=0 to 20; that is, the Θ[0,t]=t/10*2π from t=0 to t=20. As I apply this force, (some of) the other oscillators will move too as the first oscillator will exert a force on them through the spring which connects the oscillators.
• Finally, after time t=20, I no longer apply an external force on any of the oscillators and just allow it to oscillate based on the relations found in the differential equations.

If you're interested to see all the function definitions and the constants, they can be found below:

(* providing values for all our variables. these variables describe one particular set of coupled oscillators *)
M=0.1;
L=0.2;
β=1;
μ=1;
k=100;
η=0.3;
γ=0.2;
α=1;
g=9.81
(* defining the torque τ. funcf and funcg are only used to make the expression for τ simpler.*)
funcg=(1-η/(2(1-Cos[#1-#2])+(γ/α)^2)^0.5)&;
funcf=#+Abs[#]&;
τ=Sin[#1-#2]*k*L^2*funcf[funcg[#1,#2]]&;


Edit: Since I was solving a physics question, I forgot to state that what I meant by g in the equations was the gravitational acceleration on earth 9.81 and nothing related to the function g which I was using as an intermediate function to make my expression of τ less complicated. I've modified this by writing the function g as funcg.

After some work, I realized that the method suggested below of using nodes helped to simplify the code a great deal. Thanks, @Spawn1701D.

• @Nasser sorry for rushing off the question. I hope it's more understandable now with the edits, especially the part about the initial condition. Mar 11 '13 at 14:07
• If it isn't, please tell me; i'm quite sure a lot of people were confused with the initial phrasing of the question. Mar 11 '13 at 14:07
• Can you explain what you mean by the while loop? I assume that you create the grid with it but this can be easily done with the Table command for instance. Mar 11 '13 at 18:27
• g is a two-variable function on your system looks like a constant can you correct that? Mar 11 '13 at 21:17
• @VincentTjeng: I'm still not sure whether I understand your initial conditions correct: If what you are saying is that the angle of the boundary node is a predefined "forced" function, then I think you'd best split your total time interval into two parts: for the first part solve the system of n-1 oscillators and Θ[0,t]=t/10*2π fixed (Θ[0,t] then isn't a dependent variable anymore). Then use the final states of that solution as initial conditions of another call to NDSolve for the full system with n oscillators for 20<t<tmax. Does that make sense? Mar 12 '13 at 22:43

Here is something that works with the definitions you gave and does what I understand you try to do:

n = 10;

fullSystem = Join[
{
β*M*L^2*D[Θ[i][t],{t,2}] == M*g*μ*L*(-Sin[Θ[i][t]]) + τ[Θ[i+1][t],Θ[i][t]]
} /. i -> 0,
Table[
β*M*L^2*D[Θ[i][t],{t,2}]
==
M*g*μ*L*(-Sin[Θ[i][t]]) + τ[Θ[i-1][t],Θ[i][t]] + τ[Θ[i+1][t],Θ[i][t]],
{i,n-1}
],
{
β*M*L^2*D[Θ[i][t],{t,2}] == M*g*μ*L*(-Sin[Θ[i][t]]) + τ[Θ[i-1][t],Θ[i][t]]
} /. i -> n
]

Θfixed = Function[{t}, t/10*2*Pi]

sysFirstPart = Rest[fullSystem] /. Θ[_?(# == 0 &)] -> Θfixed

initFirstPart = Table[{Θ[i] == 0,Derivative[Θ[i]] == 0}, {i, 1, n}]

solFirstPart = NDSolve[
Flatten[{sysFirstPart, initFirstPart}],
Table[Θ[i], {i,n}], {t, 0, 20}, MaxSteps -> 100000
]

Plot[Evaluate[Table[Θ[i][t],{i,0,n}] /. solFirstPart], {t, 0, 20}]

initSecondPart = Join[
{Θ == Θfixed, Derivative[Θ] == Derivative[Θfixed]},
Table[Θ[i],{i,n}]
==
Table[Θ[i] /. solFirstPart[],{i,n}]
],
Table[Derivative[Θ[i]],{i,n}]
==
Table[Derivative[Θ[i]] /. solFirstPart[],{i,n}]
]
]

solSecondPart = NDSolve[
Flatten[{fullSystem, initSecondPart}],
Table[Θ[i],{i,0,n}], {t,20,30}, MaxSteps -> 100000
]

Plot[Evaluate[Table[Θ[i][t],{i,0,n}] /. solSecondPart],{t,20,30}]


I have changed the naming of the variables slightly from Θ[i,t] to Θ[i][t] which makes some of the manipulations somewhat simpler, you might have a look at this answer why I think it is a good idea to be somewhat pedantic about how to name things. Other than that I'm solving the system in two parts: for the time interval 0<=t<=20 I'm solving without the first (i==0) oscillator and insert the forced angle for which I have defined the function Θfixed. Then I determine the values of Θ[i][t] and their derivatives at the end of that time interval from that solution solFirstPart and use these as the initial conditions initSecondPart for the remaining time span for which I chose 20<=t<=30. For that second time interval I solve the full system, now including the equation and angle-position of the first oscillator. You could, if necessary, use any of the techniques mentioned in this question to "concatenate" the two solutions into one. It should also be no problem to combine this solution with the more elegant way to define the system with the nodes function that was given in another answer by Spawn1701D.

With version 9 it might be possible to solve the above in one call to NDSolve by making use of the new WhenEvent: I think you would need to reformulate your problem at least for the first oscillator so it uses only first derivatives of angle Θ[t] and angular velocity vΘ[t]. Then it should be possible to switch between the two parts in both the differential equations for Θ[t] and for vΘ[t] with a WhenEvent. I haven't tried this but if the OP or anyone else is interested and finds time to try that I'd be very interested if NDSolve will accept and solve such a system.

Note: For most applications I would expect such a one-stage solution using WhenEvent to be more complicated and probably less performant than the two-stage solution shown, but there might be cases where it is to be favoured, e.g. when one wants to couple such a system to other equations...

• Hi @Albert, thank you for the answer. I'll try to see whether NDSolve accepts a system as you suggested. First, though, could you help me understand how the function Θ[i][t] allows i to be a parameter and t to be a variable? For my code, I ended up using a fudge where I named the functions Θi[t] instead, which is way less elegant than your solution. Mar 13 '13 at 15:58
• @VincentTjeng: there is no magic in using expressions like Θ[i] as variable names for NDSolve, which only serve as unique labels. NDSolve will accept certain expressions besides pure Mathematica symbols as well. You just must ensure such expressions don't evaluate (that is, there are no definitions for them). As you have found, that can makes things a lot more elegant. You can even nest that pattern: NDSolve will happily accept expression like x or x["body1"] as variable names. Other things it accepts are e.g. Subscript[x,1], Subscript["x","1"] or even just "x"... Mar 13 '13 at 17:36
• thanks, that opens the possibilities for my naming of variables a lot. Mar 15 '13 at 1:17
• why does the code in this bit require one to write it as Θ[_?(# == 0 &)], rather than simply writing Θ? Mar 19 '13 at 7:44
• @VincentTjeng: It isn't necessary in this specific case, it's more a general measure of precaution: pattern matching with numeric quantities has some pitfalls, c.f. MatchQ[Θ // N, Θ], and my complicated pattern here is a leftover from an intermediate state with a N somewhere that I did remove lateron. Also depending on how that 0. was generated it might be something slightly differing from an "exact" 0., which Θ[_?(# == 0 &)] would also take care of (unlike e.g. Θ[0|0.]). Mar 19 '13 at 11:31

You can do something like the following: First define the DE in every node

node = β*M*L^2*D[Θ[i][t],{t,2}]==M*g*μ*L*(-Sin[Θ[i][t]])+τ[Θ[i+1][t],Θ[i][t]];
node[n] = β*M*L^2*D[Θ[i][t],{t,2}]==M*g*μ*L*(-Sin[Θ[i][t]])+τ[Θ[i-1][t],Θ[i][t]];
node[i_Integer/;1$<$i$<$n] = β*M*L^2*D[Θ[i][t],{t,2}]==M*g*μ*L*(-Sin[Θ[i][t]])+τ[Θ[i-1][t],Θ[i][t]]+τ[Θ[i+1][t],Θ[i][t]];


Then create the system of DE

system = Table[node[i],{i,1,n}]


and then solve it

sol=NDSolve[Join[system, InitialConditions],Through[Array[Θ, 20][t]],{t,0,20}]

• thanks for your answer; i'll need to read through the code first to understand it! For the while loop, I set up multiple lists containing the data, let me post it up soon. Mar 12 '13 at 0:52
• @VincentTjeng also make more clear the action of the function g. Mar 12 '13 at 7:08
• are the two g symbols clearer now? One g was intended to be used as the acceleration due to the earth's gravity, while the other g was intended just as a function that makes the expression for τ less complicated. I've renamed the function as funcg. Mar 12 '13 at 15:42
• yes now almost everything is clear! ;) Mar 12 '13 at 17:49