Final state behaviour depending on parameters in coupled ODEs

Question

I have the following set of coupled differential equations subject to some initial conditions which I can solve using NDSolve (after assigning numerical values for remaining arbitrary coefficients):

{
    Mt x''[t] == 
     - k1 x[t] - k2 (x[t] - y[t]) + 
      mc rcm Sin[ϕ[t]] ϕ''[t] - γ x'[t],

    Mt y''[t] == 
     - k1 y[t] - k2 (y[t] - x[t]) + 
      mc rcm Sin[ψ[t]] ψ''[t] - γ y'[t],

    Inertia ϕ''[t] + 
     mc g rcm Sin[ ϕ[t]] + 
      γm ϕ'[t] ((ϕ[t]/θ0)^2 - 1) + 
       mc x''[t] rcm Cos[ϕ[t]] == 0 ,

    Inertia ψ''[t] + 
     mc g rcm Sin[ ψ[t]] + 
      γm ψ'[t] ((ψ[t]/θ0)^2 - 1) + 
       mc y''[t] rcm Cos[ψ[t]] == 0
 }

subject to the following initial conditions

{
    x[1] == 1, 
    x'[1] == 1, 
    y[1] == 1, 
    y'[1] == 1, 
    ϕ[1] == 1, 
    ϕ'[1] == 1, 
    ψ[1] == 1, 
    ψ'[1] == 1
 }

The system can be solved using NDSolve. Now I have two problems.

1. Instead of a single ϕ and ψ, I want to write down and find the solution for general ϕi and ψi with remaining parameters in the equations unchanged, and where i goes from 1 to N (i.e. instead of 2 coupled ODEs involving ϕi and ψi, we have 2N coupled ODEs each involving a specific ϕi and ψi). Here N is input by hand before finding the solution. The initial conditions can also be defined accordingly. How do I write and solve for the same?

2. After this is solved, starting from the same initial conditions, I want to vary the values of k1, k2 and rcm with some step size each time; and want to get the final values of ϕi and ψi after some finite time. How do I do that?

Any suggestions will be greatly helpful.

Answer 1

I suggest you use matrix forms to generalize your equations. since I do not have general equations for xn, yn,ϕn,ψi, can only generalize what you have given to me.

I start with the following vectors:

R = {x[t], y[t], \[Phi][t], \[Psi][t]};
Q = {0, 0, \[Phi][t], \[Psi][t]};
zero = {0, 0, 0, 0};

Then I define the following matrices:

A = {{mt, 0, -mc rcm Sin[\[Phi][t]], 0}, {0, mt, 
0, -mc rcm Sin[\[Psi][t]]}, {-mc rcm Cos[\[Phi][t]], 0, Inertia, 
0}, {0, -mc rcm Cos[\[Psi][t]], 0, Inertia}};
B = {{\[Gamma], 0, 0, 0}, {0, \[Gamma], 0, 0}, {0, 
0, \[Gamma]m ((\[Phi][t]/\[Theta]0)^2 - 1), 0}, {0, 0, 
0, \[Gamma]m ((\[Psi][t]/\[Theta]0)^2 - 1)}};
Cc = {{k1 + k2, -k2, 0, 0}, {-k2, k1 + k2, 0, 0}, {0, 0, 0, 0}, {0, 0,
 0, 0}};
Dd = mc rcm Sin[Q];

The final equations can be derived and solved as follows:

eq=A.D[R, {t, 2}] + B.D[R, t] + Cc.R + Dd == zero;
NDSolve[{eq, conditions}, R, {t, tmin, tmax}]

you can see that the elements of following vector

G = A.D[R, {t, 2}] + B.D[R, t] + Cc.R + Dd

would yield the the equations correctly, i.e.

G[[1]]==0

gives the correct equation for x[t], and so on. As next step, generalizing these matrixes would be easy since you only need to add terms regarding new ϕs and ψs.