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 am trying to solve system of differential equations below

$$x''(t)=-2.25cos(1.5t)-x(t)-x'(t)$$

$$y''(t)=-1.125cos(1.5t)-4y(t)-y'(t)$$

$$z'(t)=\begin{cases}x(t)-(1+y(t)), & \text{if } x>(1+y) \\ \\ 0, & \text{if } |x|<(1+y) \\ \\ x(t)+(1+y(t)), & \text{if } x<-(1+y) \end{cases} $$

I can solve the equations for $x$ and $y$, but I have no idea how to implement WhenEvent Method to solve for $z$. Notice that there are 3 conditions for $z$.

EDIT 11/14:

I tried code below

sol1 = First@NDSolve[{x''[t] == -2.25 Cos[1.5 t] - x[t] - x'[t],
 x[0] == 0, x'[0] == 0}, x, {t, 0, 30}]
sol2 = First@NDSolve[{y''[t] == -1.125 Cos[1.5 t] - 4 y[t] - y'[t],
 y[0] == 0, y'[0] == 0}, y, {t, 0, 30}]

sol3 = NDSolve[z'[t] == a[t], z[0] = 0,
 WhenEvent[Evaluate[x[t] /. sol1] > 1 + Evaluate[y[t] /. sol2],
  a[t] -> Evaluate[x[t] /. sol1 - (1 + Evaluate[y[t] /. sol2])]],
 WhenEvent[Abs[Evaluate[x[t] /. sol1]] < 1 + Evaluate[y[t] /. sol2],
  a[t] -> 0],
 WhenEvent[Evaluate[x[t] /. sol1] < -(1 + Evaluate[y[t] /. sol2]),
  a[t] -> Evaluate[x[t] /. sol1] + (1 + Evaluate[y[t] /. sol2])], 
z, {t, 0, 30}]

Plot[Evaluate[x[t] /. sol1], {t, 0, 10}, PlotRange -> All]
Plot[Evaluate[y[t] /. sol2], {t, 0, 10}, PlotRange -> All]
Plot[Evaluate[z[t] /. sol3], {t, 0, 10}, PlotRange -> All]

But it returns many errors I don't understand.
What is the proper way to solve this problem?

share|improve this question

1 Answer 1

up vote 2 down vote accepted

You can use Piecewise function in the definition of DE. So for system at hand the NDSolve command would be

xyz = First@
  NDSolve[{x''[t] == -2.25 Cos[1.5 t] - x[t] - x'[t], x[0] == 0, 
    x'[0] == 0, y''[t] == -1.125 Cos[1.5 t] - 4 y[t] - y'[t], 
    y[0] == 0, y'[0] == 0, 
    z'[t] == 
     Piecewise[{{x[t] - (1 + y[t]), x[t] > 1 + y[t]}, {0, 
        Abs[x[t]] <= 1 + y[t]}, {x[t] + (1 + y[t]), 
        x[t] < -1 - y[t]}}], z[0] == 0}, {x, y, z}, {t, 0, 30}];

Then the plot of the solution is

Plot[{x[t], y[t], z[t]} /. xyz // Evaluate, {t, 0, 30}, 
 PlotLegends -> {"x(t)", "y(t)", "z(t)"}, AxesLabel -> {"t", None}]

enter image description here

share|improve this answer
    
Very nice! Thank you mmal! I didn't know you can use piecewise inside NDSolve! And also, if I may ask, what is the function of "First@" before NDSolve? I often see that, but couldn't find the tutorial. –  Dadan Ari Wibowo Nov 14 '13 at 15:27
    
You're welcome. I have used First since NDSolve returns list of solutions (in general there may be more than one solution), I'm taking first element of that list (which in this particular case has exactly one element). Just compare output of NDSolve[...] with First[NDSolve[...]] and how to use them, etc. –  mmal Nov 15 '13 at 9:45

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.