Message Boards Message Boards

0
|
248 Views
|
3 Replies
|
2 Total Likes
View groups...
Share
Share this post:

Error 1/0 encountered while plotting solution of differential equations

Posted 13 days ago

Hi everyone. I've started to use Wolfram Mathematica recently and I was trying to solve a system of coupled differential equations. The second of these equations presents a function dependent by the "t " variable and defined in a specific range of values, since it presents points of singularity. Unfortunately, even by defining all the intervals where such function is defined, it seems like the solver keeps considering these points belonging to the interval, giving the message "ERROR: infinite expression 1/0 encountered".

I'll write in the following lines what I've put in the program:

m1 = 50; (*kg*)
m2 = 8; (*kg*)
k1 = 4000; (*N/m*)
k2 = 180000; (*N/m*)
c = 150; (*Ns/m*)
v = 2.78; (*m/s*)
z = 2.71828;
H0 = 0.09; (*m*)
L0 = 0.27; (*m*)


ClearAll[h];
h[t_] := 
   Piecewise[{{0, 
      t <= 1}, {Unevaluated[
       z*H0*Exp[-1/(1 - (2*v*(t - 1)/L0 - 1)^2)]], 
      1 < t < 1 + L0/v}, {0, t >= 1 + L0/v}}]; // FullSimplify


sol = NDSolve[{m1  x''[t] + c  (x'[t] - y'[t]) + k1  (x[t] - y[t]) == 
     0, m2  y''[t] + c  (y'[t] - x'[t]) + k1  (y[t] - x[t]) + 
      k2  (y[t] - h[t]) == 0, x[0] == 0, x'[0] == 0, y[0] == 0, 
    y'[0] == 0}, {x[t], y[t]}, {t, 0, 4.5}];

Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 4.5}, 
 PlotLabels -> {"x1(t)", "x2(t)"}, 
 AxesLabel -> {"time, s", "displacement, m"}, 
 PlotLegends -> "Expressions", PlotRange -> All, 
 Exclusions -> "Singularities"]

I've already tried to check in another discussions about this topic, but the problem sill persists. So the final question is, how can I solve the error here presented? Thank you!

Attachment

Attachments:
3 Replies

You're welcome.

POSTED BY: Michael Rogers

The problem is with NDSolve[], not with plotting per se. I think it's probably a bug in handling the discontinuity in Piecewise. See if this workaround is satisfactory:

sol = NDSolve[{m1   x''[t] + c   (x'[t] - y'[t]) + 
      k1   (x[t] - y[t]) == 0, 
    m2   y''[t] + c   (y'[t] - x'[t]) + k1   (y[t] - x[t]) + 
      k2   (y[t] - h[t]) == 0, x[0] == 0, x'[0] == 0, y[0] == 0, 
    y'[0] == 0} // Simplify`PWToUnitStep, {x[t], y[t]}, {t, 0, 4.5}]

Alternatively, we can manually process the discontinuities:

sol = NDSolve[{m1   x''[t] + c   (x'[t] - y'[t]) + 
     k1   (x[t] - y[t]) == 0, 
   m2   y''[t] + c   (y'[t] - x'[t]) + k1   (y[t] - x[t]) + 
     k2   (y[t] - h[t]) == 0, x[0] == 0, x'[0] == 0, y[0] == 0, 
   y'[0] == 0
   , WhenEvent[t == 1, "RestartIntegration"]
   , WhenEvent[t == 1 + L0/v, "RestartIntegration"]
   }, {x[t], y[t]}, {t, 0, 4.5}, 
  Method -> "DiscontinuityProcessing" -> None]
POSTED BY: Michael Rogers

Yes, it finally works. I don't know how to thank you properly! Thank you.

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