Message Boards Message Boards

GROUPS:

Solve a second order ODE modeling a mass-spring system?

Posted 2 months ago
437 Views
|
7 Replies
|
5 Total Likes
|

I am new to Mathematica and still learning the language. I have question involving a mass -spring system and need to solve a second order ODE for multiple points at multiple times but can't write it properly. Please help the question and code written is attached.

g = 10

x = Range[-.5, .5, 1/13]

DSolve[p''[x,y[t]] == -100 (Norm[{p[x, y[t]]} - {p[x - 1, 
         y[t]]}])*(UnitVector[{p[x, y[t]]} - {p[x - 1, y[t]]}]) - 
   100 (Norm[{p[x, y[t]]} - {p[x + 1, y[t]]}])*(UnitVector[{p[x, 
         y[t]]} - {p[x + 1, y[t]]}]) - 
   10 (p'[x, y[t]] - p'[x - 1, y[t]]) - 
   10 (p'[x, y[t]] - p'[x + 1, y[t]]), p[x, y[t]], y[t]]
7 Replies

You should post your code using the code sample icon (the first one).

Posted 2 months ago

Perhaps starting with a slightly simplified version of your problem with a few small errors corrected and try to figure out why it is not able to understand a couple of the things you have done

g = 10;

system = Table[
  p''[x, y[t]] ==
   -100 Norm[{p[x,y[t]]}-{p[x-1/3,y[t]]}]*UnitVector[{p[x,y[t]]}-{p[x-1/3,y[t]]}] -
    100 Norm[{p[x,y[t]]}-{p[x+1/3,y[t]]}]*UnitVector[{p[x,y[t]]}-{p[x+1/3,y[t]]}] -
    10 (p'[x,y[t]]-p'[x-1/3,y[t]])-10 (p'[x,y[t]] - p'[x+1/3,y[t]])+g,
    {x, -1/2 + 1/3, 1/2 - 1/3, 1/3}
]

functions = Table[p[x, y[t]], {x, -1/2, 1/2, 1/3}]

DSolve[system, functions, y[t]]

The Abs and UnitVector which remain after it has finished calculating are an indication that it doesn't know how to apply these to your p and y functions. Perhaps you can find another way of expressing your problem.

If you can resolve those problems then you should be able to replace 1/3 with 1/13 and solve your full problem.

Posted 2 months ago

Thanks for the simplification. I understood the need for the simplification and can represent the acceleration and velocity as pure y''[t] and y'[t] due to symmetry but the position is needed to find the extension in the spring. I cant figure out another way to solve the extension in the spring.

Posted 2 months ago

[EDIT]Made it simpler after a few assumptions. Mathematica has been running for 5 minutes without producing an answer. HELP

g = 10;
ic1 = {y[x, 0] == 0; D[y[x, 0], {t, 1}] == 0}
ic2 = {y[x, 0] ==  0; y'[x, 0] == -100 x^2 + 25}
ic3 = {y[x, 0] ==  
Which[x < 0, y[x, 0] = 1/2 x + 1/4, x > 0, 
y[x, 0] = -1/2 x + 1/4 ]; y'[x, 0] == 0}
bc = {y[-1/2, t] == 0, D[y[-1/2, t] == 0, {t, 1}], 
D[y[-1/2, t] == 0, {t, 2}], y[1/2, t] == 0; 
D[y[1/2, t] == 0, {t, 1}] 0, D[y[1/2, t] == 0, {t, 2}]}
system = Table[
  D[y[x, t], {t, 2}] == -100 (y[x, t] - y[x - 1/13, t]) - 
    100 (y[x, t] - y[x + 1/13, t]) - 
    10 (D[y[x, t], {t, 1}] - D[y[x - 1/13, t], {t, 1}]) - 
    10 (D[y[x, t], {t, 1}] - D[y[x + 1/13, t], {t, 1}]) + 
    g, {x, -1/2 + 1/13, 1/2 - 1/13, 1/13}]

functions = [y[x, t], {x, -1/2, 1/2, 1/13}]

DSolve[{system, bc, ic1}, functions, t]

I will show the working code and the solution to task 2

t0 = 9; P = 14; lr = 1; g = 9.8; m = 1; ks = 100; kd = 10;
eqns = {Table[
    m*x[i]''[t] == -ks*(Sqrt[(x[i][t] - x[i - 1][t])^2 + (y[i][t] - 
              y[i - 1][t])^2] - lr)*(x[i][t] - x[i - 1][t])/
        Sqrt[(x[i][t] - x[i - 1][t])^2 + (y[i][t] - y[i - 1][t])^2] - 
      ks*(Sqrt[(x[i][t] - x[i + 1][t])^2 + (y[i][t] - 
              y[i + 1][t])^2] - lr)*(x[i][t] - x[i + 1][t])/
        Sqrt[(x[i][t] - x[i + 1][t])^2 + (y[i][t] - y[i + 1][t])^2] - 
      kd*(x[i]'[t] - x[i - 1]'[t]) - kd*(x[i]'[t] - x[i + 1]'[t]), {i,
      2, P - 1}], 
   Table[m*y[i]''[
       t] == -ks*(Sqrt[(x[i][t] - x[i - 1][t])^2 + (y[i][t] - 
              y[i - 1][t])^2] - lr)*(y[i][t] - y[i - 1][t])/
        Sqrt[(x[i][t] - x[i - 1][t])^2 + (y[i][t] - y[i - 1][t])^2] - 
      ks*(Sqrt[(x[i][t] - x[i + 1][t])^2 + (y[i][t] - 
              y[i + 1][t])^2] - lr)*(y[i][t] - y[i + 1][t])/
        Sqrt[(x[i][t] - x[i + 1][t])^2 + (y[i][t] - y[i + 1][t])^2] - 
      kd*(y[i]'[t] - y[i - 1]'[t]) - kd*(y[i]'[t] - y[i + 1]'[t]) - 
      m*g, {i, 2, P - 1}]};
ic = {Table[y[i]'[0] == 0, {i, 1, P}], Table[y[i][0] == 0, {i, 1, P}],
    Table[x[i][0] == lr*(i - 1), {i, 1, P}], 
   Table[x[i]'[0] == 0, {i, 1, P}]};
bc = {x[P]''[t] == 0 , x[1]''[t] == 0, y[P]''[t] == 0, y[1]''[t] == 0};

Y = NDSolveValue[{eqns, ic, bc}, 
  Table[y[i][t], {i, 1, P}], {t, 0, t0}]; X = 
 NDSolveValue[{eqns, ic, bc}, Table[x[i][t], {i, 1, P}], {t, 0, t0}];
{Plot[Y, {t, 0, t0}, PlotRange -> All, AxesLabel -> {"t", "y"}], 
 Plot[X, {t, 0, t0}, PlotRange -> All, AxesLabel -> {"t", "x"}]}

This figure shows the coordinates of the point masses. fig1 These coordinates can be displayed on the plane.

Table[Graphics[{PointSize[Large], Red, 
   Point[NDSolveValue[{eqns, ic, bc}, 
     Table[{x[i][tn], y[i][tn]}, {i, 1, P}], {t, 0, tn}]]}, 
  Frame -> True, PlotLabel -> tn, 
  PlotRange -> {{-.5, 13.5}, {-9, 1}}], {tn, 0, t0, 1}]

fig2

Velocity components of point masses

VY = NDSolveValue[{eqns, ic, bc}, 
  Table[y[i]'[t], {i, 1, P}], {t, 0, t0}]; VX = 
 NDSolveValue[{eqns, ic, bc}, Table[x[i]'[t], {i, 1, P}], {t, 0, t0}];
{Plot[VY, {t, 0, t0}, PlotRange -> All, AxesLabel -> {"t", "v"}], 
 Plot[VX, {t, 0, t0}, PlotRange -> All, AxesLabel -> {"t", "u"}]}

fig3

Posted 2 months ago

Thank for this amazing code. You are a God-sent.

Just one technical question. The boundary conditions are zero acceleration at the endpoints. Would Zero displacement work?

You can use any boundary conditions. This code is adapted for your task. However, you can use the code to solve a problem with one or two free ends. Therefore, accelerations are used as boundary conditions. An example of solving a problem with one free end fig4

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract