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

Question

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.

Answer 1

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

node[1] = β*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}]

Answer 2

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] == 0,Derivative[1][Θ[i]][0] == 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[
  {Θ[0][20] == Θfixed[20], Derivative[1][Θ[0]][20] == Derivative[1][Θfixed][20]},
  Thread[
      Table[Θ[i][20],{i,n}] 
      == 
      Table[Θ[i][20] /. solFirstPart[[1]],{i,n}]
  ],
  Thread[
      Table[Derivative[1][Θ[i]][20],{i,n}] 
      == 
      Table[Derivative[1][Θ[i]][20] /. solFirstPart[[1]],{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 Θ[0][t] and angular velocity vΘ[0][t]. Then it should be possible to switch between the two parts in both the differential equations for Θ[0][t] and for vΘ[0][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...