Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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?

share|improve this question
1  
Have you tried something? Can you post your (possibly failing) code? –  belisarius Oct 25 '13 at 15:57
1  
Check out the documentation for Condition –  Timothy Wofford Oct 25 '13 at 16:39
add comment

3 Answers

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

Mathematica graphics

share|improve this answer
    
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
add comment

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}

enter image description here

share|improve this answer
2  
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
add comment

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

Harmonic oscillator with various compliant stoppers

share|improve this answer
    
Thank you for your suggestion, the Piecewise function works well. –  user10204 Nov 15 '13 at 14:58
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.