# Find numerical solution to this system of DE

I am trying to solve this system

$$\left( \begin{array}{ccccc} 2 k & -k & 0 & 0 & 0 \\ -k & 2 k & -k & 0 & 0 \\ 0 & -k & 2 k & -k & 0 \\ 0 & 0 & -k & 2 k & -k \\ 0 & 0 & 0 & -k & 2 k \\ \end{array} \right) \left( \begin{array}{c} y_1(t) \\ y_2(t) \\ y_3(t) \\ y_4(t) \\ y_5(t) \\ \end{array} \right)+\left( \begin{array}{ccccc} m & 0 & 0 & 0 & 0 \\ 0 & m & 0 & 0 & 0 \\ 0 & 0 & m & 0 & 0 \\ 0 & 0 & 0 & m & 0 \\ 0 & 0 & 0 & 0 & m \\ \end{array} \right) \left( \begin{array}{c} y_1''(t) \\ y_2''(t) \\ y_3''(t) \\ y_4''(t) \\ y_5''(t) \\ \end{array} \right)=\left( \begin{array}{c} F(t) \\ 0 \\ F(t) \\ 0 \\ F(t) \\ \end{array} \right)$$ but the following code doesn't seem to work..

n = 6;

M = m IdentityMatrix[n - 1];

K = Table[0, {j, 1, n - 1}, {i, 1, n - 1}];

For[j = 1, j <= n - 1, j++,
For[i = 1, i <= n - 1,
i++, {If[j == i, K[[j, i]] = 2 k, Nothing],
If[i == j + 1, K[[j, i]] = -k, Nothing],
If[i == j - 1, K[[j, i]] = -k, Nothing]}]]

yy[t] = Table[Subscript[y, i][t], {i, 1, n - 1}];

FF[t] = Table[
If[(i/n) == 1/6 || i == n/2 || (i/n) == 5/6, F[t], 0], {i, 1,
n - 1}];

F[t_] = Piecewise[{{Subscript[F, 0] t/(T/2),
0 <= t <= T/2}, {-Subscript[F, 0],
T/2 <= t <= T}}];

m = 1;
k = 1;
Subscript[F, 0] = 1;
T = 1;

NDSolve[{M.D[yy[t], {t, 2}] + K.yy[t] == F[t], yy[t] == 0 /. t -> 0,
D[yy[t], t] == 0 /. t -> 0}, yy[t], {t, 10 T, 20 T}]

• You can do away with your For loops and have much more readable code (though your layout is good) by using DiagonalMatrix: reference.wolfram.com/language/ref/DiagonalMatrix.html Commented May 16, 2016 at 20:08
• @Quantum_Oli Wow, I wasn't aware of that option! Thank you for that! Commented May 16, 2016 at 21:00

Mathematica is checking to make sure that you have enough equations and unknowns, but the way that you've written them out, it only thinks that you have three equations. This is because you've set them up as {5-component vector} == 0, rather than {5-component vector} == {5-component vector}.

There are a couple of ways to fix this. One is to use Thread:

soln = NDSolve[{Thread[M.D[yy[t], {t, 2}] + K.yy[t] == FF[t]], Thread[yy[t] == 0] /. t -> 0, Thread[D[yy[t], t] == 0] /. t -> 0}, yy[t], {t, 10 T, 20 T}]


(N.B.: I have replaced F[t] with FF[t] here, which I suspect is what you really want given your original equations at the top.) This evaluates with no errors, giving the following result:

Plot[Evaluate[{Subscript[y, 1][t], Subscript[y, 2][t],  Subscript[y, 3][t], Subscript[y, 4][t], Subscript[y, 5][t]} /. First[soln]], {t, 10, 20}]


(Note that $y_2 = y_4$ and $y_1 = y_5$ given your initial conditions and equations, which is why only three graphs are visible in the plot below.)

Alternately, you can write out the equations with vectors on both sides:

soln = NDSolve[{M.D[yy[t], {t, 2}] + K.yy[t] == FF[t],
yy[t] == {0, 0, 0, 0, 0} /. t -> 0,
D[yy[t], t] == {0, 0, 0, 0, 0} /. t -> 0}, yy[t], {t, 10 T, 20 T}]


This yields the same result.

• I am not sure why but this evaluates with an error in both cases saying NDSolve::depdole: "The differential order of a dependent variable in (_all_derivatives__) exceeds the highest order that appears in the differential equations. " Commented May 16, 2016 at 20:59
• Hmm, I'm not sure what to tell you; both methods work fine for me with a fresh kernel. (Apologies if this is obvious, but have you tried quitting the kernel and restarting it? It often clears up lingering definitions of variables that are getting in the way.) I'm using v10.2 on Mac OS, if it matters. Commented May 16, 2016 at 21:42
• I tried quitting kernel and restarting also. Also restarted my computer in case that would make any difference. I am using version 10.3.1.0 (student version) on a microsoft windows 10 (64 bit). Well I would love to accept your answer yet it doesn't work. I apologize for that. Here is a screenshot of my code: www4.slikomat.com/13/0517/5bx-Captur.png Commented May 17, 2016 at 9:56
• @skrat: I think you've run into a known bug in v10.3, arising when discontinuous functions are fed into NDSolve. See this answer for a work-around. Commented May 17, 2016 at 13:16
• Yep, you are right. It works now! Commented May 17, 2016 at 15:37