# Differential equation of mass spring system with compliant stoppers

I would like to solve numerically the differential equation for the displacement x[t] of a mass m-spring k system with compliant stoppers.

The equation should be something like:

m x"[t] == -k x[t] -F[t]


F[t] should be defined so that:

-it is 0 when Abs[x] is below a certain value x0

-it is equal to - k2 (x[t] - x0) when Abs[x] is above x0.

Which Mathematica function should I use to define F[t] so that NDSolve can solve the corresponding equation?

-
Have you tried something? Can you post your (possibly failing) code? –  belisarius Oct 25 '13 at 15:57
Check out the documentation for Condition –  Timothy Wofford Oct 25 '13 at 16:39

k = 1;
m = 1;
x0 = 1/2;
f1[k2_][t_, x_] := -k/m x;
f2[k2_][t_, x_] := -k/m x + k2/m (x0 - x);
f[k2_][t_, x_] := If[Abs@x <= x0, f1[k2][t, x], f2[k2][t, x]];
fb[k2_][t_?NumberQ, x_] := f[k2][t, x];
bsol = NDSolveValue[{x''[t] == fb[#][t, x[t]], x[0] == 0, x'[0] == 1}, x, {t, 0, 10}] & /@ Range@5
ListLinePlot[bsol, Mesh -> All, PlotRange -> All, GridLines -> {{}, {-x0, x0}}]


-
It seems like on v9, the WhenEvent option could be used in this case to toggle the forcing function. –  Timothy Wofford Oct 25 '13 at 19:31
@TimothyWofford Yep, but WhenEvent has some subtleties that I don't want to confront today :) –  belisarius Oct 25 '13 at 19:36

You can adapt the code based on the function you have. Free vibration of SDOF systems (undamped) draws displacement velocity and acceleration with respect to time can be calculated as in below code. Some demos can be found at Wolfram demonstrations site. demo1 demo2 demo3

m = 2;
k = 8;
w = Sqrt[k/m];

u = DSolve[{y''[x] + (w^2)*y[x] == 0, y[0] == 3, y'[0] == 5}, y[x], x][[1, 1, 2]] // Simplify


3 Cos[2 x] + 5 Cos[x] Sin[x]

v = D[u, x]


5 Cos[x]^2 - 5 Sin[x]^2 - 6 Sin[2 x]

a = D[v, x]


-12 Cos[2 x] - 20 Cos[x] Sin[x]

p1 = Plot[u, {x, 0, 2}, AxesLabel -> {"time in secs", "displacement"},
PlotStyle -> {Thickness[0.01]}];
p2 = Plot[v, {x, 0, 2}, AxesLabel -> {"time in secs", "velocity"},
PlotStyle -> {Thickness[0.01]}];
p3 = Plot[a, {x, 0, 2}, AxesLabel -> {"time in secs", "acceleration"},
PlotStyle -> {Thickness[0.01]}];
{p1, p2, p3}


-
I guess you missed the requirement with respect to the constraints on F? –  Sjoerd C. de Vries Oct 25 '13 at 17:39
Yes you are right. The Question was not so clear. So, I gave a sample code to start with toying. –  s.s.o Oct 25 '13 at 17:48

Belisarius has a nice solution using If which I am more or less copying here for comparison. I've removed t from the forcing function because it doesn't actually depend on time explicitly. Note, both of us have assumed that you meant F[t]==+k2(x-x0) when Abs@x > 0 because otherwise you get exponential growth.

k = 1; m = 1; x0 = 1/2;
F0[k2_][x_] := If[Abs@x > x0, k2 (x - x0), 0]


In the comments I suggested Condition

F1[k2_][x_] /; Abs@x <= x0 := 0
F1[k2_][x_] /; Abs@x > x0 := k2 (x - x0)


Another option is Piecewise

F2[k2_][x_] := Piecewise[{{k2 (x - x0), Abs@x > x0}}, 0]


And for those of us still running v8,

bsol = x[t] /. NDSolve[{
m x''[t] == -k x[t] - F2[#][x[t]],
x[0] == 0, x'[0] == 1}, x, {t, 0, 10}] & /@ Range@5;
Plot[bsol, {t, 0, 10}]


-
Thank you for your suggestion, the Piecewise function works well. –  user10204 Nov 15 '13 at 14:58