Message Boards

WOLFRAM COMMUNITY

7900 Views

7 Replies

5 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Control Systems Engineering Mathematics Physics Equation Solving Wolfram Language Modeling

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

Areeb Qureshi

Posted 6 years ago

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]]

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]]

POSTED BY: Areeb Qureshi

7 Replies

Sort By:

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 6 years ago

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

POSTED BY: Frank Kampas

Bill Nelson

Posted 6 years 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 BY: Bill Nelson

Areeb Qureshi

Posted 6 years 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 BY: Areeb Qureshi

Areeb Qureshi

Posted 6 years 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]

[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]

POSTED BY: Areeb Qureshi

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 6 years ago

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[ mx[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[my[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]) - mg, {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. 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}] 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"}]}

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 BY: Alexander Trounev

Areeb Qureshi

Posted 6 years 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?

POSTED BY: Areeb Qureshi

Alexander Trounev

Alexander Trounev, A&E Trounev IT Consulting

Posted 6 years ago

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

POSTED BY: Alexander Trounev

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback