当前位置：首页 > news >正文

在Mathematica中使用WhenEvent求解微分方程

news 2025/6/4 9:09:41

WhenEvent[event,action]指定当事件event触发时，方程在 NDSolve 及相关函数中执行的操作action。

模拟一个每次弹起后保持95%速度的弹跳球

NDSolve[{y''[t] == -9.81, y[0] == 5, y'[0] == 0, WhenEvent[y[t] == 0, y'[t] -> -0.95 y'[t]]}, y, {t, 0, 10}]; Plot[y[t] /. %, {t, 0, 10}]

模拟沿阶梯弹跳下落的球

c = .75;
sol = NDSolve[{y''[t] == -9.8, y[0] == 13.5, y'[0] == 5, a[0] == 13, WhenEvent[y[t] - a[t] == 0, y'[t] -> -c y'[t]], WhenEvent[Mod[t, 1], a[t] -> a[t] - 1]}, {y, a}, {t, 0, 8}, DiscreteVariables -> {a}];
Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, Filling -> {2 -> 0}]

绘制球体的动能、势能及总能量曲线图

kin[v_] := .5 v^2;
pot[y_] := 9.8 y;
energy[y_, v_] := kin[v] + pot[y];
Plot[Evaluate[{kin[y'[t]], pot[y[t]], energy[y[t], y'[t]]} /. sol], {t, 0, 8}]

在方形盒子中，建立球体碰撞侧壁时改变方向的运动模型

sol = NDSolve[{x'[t] == a[t], y'[t] == b[t], x[0] == 0, y[0] == 0, a[0] == 1, b[0] == Rationalize[Sqrt[2], .01], WhenEvent[x[t]^2 == 1, a[t] -> -a[t]], WhenEvent[y[t]^2 == 1, b[t] -> -b[t]]}, {x, y}, {t, 0, 100}, DiscreteVariables -> {a, b}];ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False]

当初速度向量为无理数时，系统不再存在周期性解

sol = NDSolve[{x'[t] == a[t], y'[t] == b[t], x[0] == 0, y[0] == 0, a[0] == 1, b[0] == Sqrt[2], WhenEvent[x[t]^2 == 1, a[t] -> -a[t]],WhenEvent[y[t]^2 == 1, b[t] -> -b[t]]}, {x, y}, {t, 0, 100}, DiscreteVariables -> {a, b}];
ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False]