Wolfram 语言与系统参考资料中心

WhenEvent

WhenEvent[event,action]

指定在 NDSolve 中，当方程中触发 event 产生的一个 action 和相关函数.

范例

打开所有单元关闭所有单元

基本范例 (1)

模拟一个在每次反弹时仍保留95%速度的反弹球：

Wolfram Language code: NDSolve[{y''[t] == -9.81, y[0] == 5, y'[0] == 0, WhenEvent[y[t] == 0, y'[t] -> -0.95y'[t]]}, y, {t, 0, 10}];

Wolfram Language code: Plot[y[t] /. %, {t, 0, 10}]

范围 (16)

事件 (6)

当 x[t] 穿过 1/2 时，输出 t 的值：

Wolfram Language code: NDSolve[{x''[t] == -t x[t], x[0] == 0, x'[0] == 1, WhenEvent[x[t] == 1 / 2, Print[t]]}, x, {t, 0, 6}];

Wolfram Language code: Plot[{x[t] /. %, 1 / 2}, {t, 0, 6}]

当 x[t] 大于 1/2 时，输出 t 的值：

Wolfram Language code: NDSolve[{x''[t] == -t x[t], x[0] == 0, x'[0] == 1, WhenEvent[x[t] > 1 / 2, Print[t]]}, x, {t, 0, 6}];

当 x[t] 小于 1/2时，输出 t 的值：

Wolfram Language code: NDSolve[{x''[t] == -t x[t], x[0] == 0, x'[0] == 1, WhenEvent[x[t] < 1 / 2, Print[t]]}, x, {t, 0, 6}];

按固定时间间隔，采样解：

Wolfram Language code:

Manipulate[{sol, samples} = 
	Reap@NDSolve[{x'[t] == Sin[Exp[t]], x[0] == 0, WhenEvent[Mod[t, Δt] == 0, Sow[{t, x[t]}]]}, x, {t, 3}];
	Plot[x[t] /. sol, {t, 0, 3}, Epilog -> {PointSize[Medium], Red, Point@@samples}], {{Δt, 1}, .01, 2}, SaveDefinitions -> True]

以列表指定多个事件，当发生任一事件便触发行动（action）：

Wolfram Language code:

NDSolve[{x''[t] + x[t] == 0, x'[0] == π, x[0] == 0, WhenEvent[{x[t] == 2, x[t] == -3}, x'[t] -> -x'[t]]}, x, {t, 0, 20}];

Wolfram Language code: Plot[{x[t] /. %, 2, -3}, {t, 0, 20}]

使用具有条件的事件检测正轴的相交点：

Wolfram Language code:

{sol, points} = Reap@NDSolve[{x'[t] == -y[t] - x[t] ^ 2, y'[t] == 2x[t] - y[t] ^ 3, x[0] == y[0] == 1, WhenEvent[y[t] == 0 && x[t] > 0, Sow[{x[t], y[t]}]]}, {x, y}, {t, 14}];

Wolfram Language code:

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 14}, Epilog -> {PointSize[Medium], Red, Point@@points}, AxesLabel -> {x, y}]

一般而言，事件可以是任何可以被检验为 True 的表达式：

Wolfram Language code:

Block[{eventtime = AbsoluteTime[] + 0.1}, 
	NDSolve[{x''[t]  + t x[t] == 0, x[0] == 1, x'[0] == 0, WhenEvent[AbsoluteTime[] > eventtime, eventtime += 0.3;Print[{AbsoluteTime[], t, x[t]}], "LocationMethod" -> "StepEnd"]}, x, {t, 0, 100}, WorkingPrecision -> 20]]

当点击按钮时，触发一个事件停止积分：

Wolfram Language code: Button["Stop", stop = True]

如果很快点击 Stop 按钮，积分会在前终止：

Wolfram Language code:

stop = False;usol = NDSolveValue[{D[u[t, x], t, t] == D[u[t, x], x, x], u[0, x] == Exp[-10 x ^ 2], Derivative[1, 0][u][0, x] == 0, u[t, -10] == u[t, 10], WhenEvent[stop, "StopIntegration"]}, u, {t, 0, 100}, {x, -10, 10}, Method -> "StiffnessSwitching"]

行为 (9)

以列表指定多个行为，按连续顺序执行：

Wolfram Language code:

NDSolve[{x''[t] == 0, x'[0] == 1, x[0] == 0, WhenEvent[{x[t] == 2, x[t] == -3}, {Print["Before: t = ", t, " x[t] = ", x[t], " x'[t] = ", x'[t]], x'[t] -> -x'[t], Print["After: t = ", t, " x[t] = ", x[t], " x'[t] = ", x'[t]]}]}, x, {t, 0, 20}];

Wolfram Language code: Plot[{x[t] /. %, 2, -3}, {t, 0, 20}]

使用 "StopIntegration" 结束一个事件的解：

Wolfram Language code: NDSolve[{x''[t] + x[t] == 0, x[0] == x'[0] == 1, WhenEvent[x[t] == 0, end = t;"StopIntegration"]}, x, {t, ∞}];

Wolfram Language code: Plot[x[t] /. %, {t, 0, end}]

在事件第一次发生后，使用 "RemoveEvent" 删除事件：

Wolfram Language code:

NDSolve[{x''[t] + x[t] == 0, x[0] == x'[0] == 1, WhenEvent[x[t] == 0, {x'[t] -> -x'[t], "RemoveEvent"}]}, x, {t, 0, 10}];

Wolfram Language code: Plot[x[t] /. %, {t, 0, 10}]

在事件第三次发生后，删除事件：

Wolfram Language code:

NDSolve[{x''[t] + x[t] == 0, x[0] == x'[0] == 1, c[0] == 0, WhenEvent[x[t] == 0, If[c[t] < 3, {x'[t], c[t]} -> {-x'[t], c[t] + 1}, "RemoveEvent"]]}, x, {t, 0, 6Pi}, DiscreteVariables -> Element[c, Integers]];

Wolfram Language code: Plot[x[t] /. %, {t, 0, 6Pi}]

在事件后，使用 "RestartIntegration" 强迫重新开始积分：

Wolfram Language code: blackbox[x_ ? NumberQ] := a x;

Wolfram Language code:

xsol = Block[{a = 1}, NDSolve[{x'[t] == blackbox[x[t]], x[0] == 1, WhenEvent[x[t] ^ 2 - 3x[t] + 1 == 0, a = -a;"RestartIntegration"]}, x, {t, 0, 5}]];

Wolfram Language code: Plot[x[t] /. xsol, {t, 0, 5}]

处理这种变化的更好的方法是使用离散状态变量 a[t]：

Wolfram Language code:

xsold = NDSolve[{x'[t]  == a[t] x[t], x[0] == 1, a[0] == 1, WhenEvent[x[t] ^ 2 - 3x[t] + 1 == 0, a[t] -> -a[t]]}, x, {t, 0, 5}, DiscreteVariables -> a];

Wolfram Language code: Plot[x[t] /. xsold, {t, 0, 5}]

规则可用于修改 DependentVariables 或 DiscreteVariables：

Wolfram Language code:

NDSolve[{y'[t] == a[t] y[t], a[0] == 5, y[0] == 1, WhenEvent[y[t] == 10, {y[t] -> -10, a[t] -> -5}]}, {y, a}, {t, 0, 1}, DiscreteVariables -> {a}  ];

Wolfram Language code: Plot[Evaluate[{y[t], a[t]} /. %], {t, 0, 1}]

使用 "CrossDiscontinuity" 管理不连续相交点：

Wolfram Language code: blackbox[x_ ? NumberQ] := If[x > 1, -1, -3];

Wolfram Language code: NDSolve[{x'[t] == blackbox[x[t]], WhenEvent[x[t] == 1, "CrossDiscontinuity"], x[0] == 2}, x, {t, 0, 2}];

Wolfram Language code: Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}]

使用 "CrossSlidingDiscontinuity" 表明可以使用菲利波夫（Filippov）延续：

Wolfram Language code: blackbox[x_ ? NumberQ] := If[x > 1, -1, 1];

Wolfram Language code: NDSolve[{x'[t] == blackbox[x[t]], x[0] == 2, WhenEvent[x[t] == 1, "CrossSlidingDiscontinuity"]}, x, {t, 0, 2}];

Wolfram Language code: Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}]

使用 DiscreteVariables 和 "DiscontinuitySignature" 有效穿越不连续点：

Wolfram Language code:

NDSolve[{x'[t] == v[t] - 2, x[0] == 2, v[0] == 1, WhenEvent[x[t] == 1, v[t] -> "DiscontinuitySignature"]}, {x, v}, {t, 0, 2}, DiscreteVariables -> {Element[v, {-1, 1}]}];

Wolfram Language code: {Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}], Plot[v[t] /. % , {t, 0, 2}]}

在 Filippov 滑动模式情况下，"DiscontinuitySignature" 为 0：

Wolfram Language code:

NDSolve[{x'[t] == -v[t], x[0] == 2, v[0] == 1, WhenEvent[x[t] == 1, v[t] -> "DiscontinuitySignature"]}, {x, v}, {t, 0, 2}, DiscreteVariables -> {Element[v, {-1, 0, 1}]}];

Wolfram Language code: {Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}], Plot[v[t] /. % , {t, 0, 2}]}

求解器 (1)

WhenEvent 可用于 NDSolve:

Wolfram Language code: ndsol = NDSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-1) ^ t]}, x, {t, 0, 3}]

Wolfram Language code: Plot[x[t] /. ndsol[[1]], {t, 0, 3}]

DSolve:

Wolfram Language code: dsol = DSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-1) ^ t]}, x, {t, 0, 3}]

Wolfram Language code: Plot[x[t] /. dsol[[1]], {t, 0, 3}]

ParametricNDSolve:

Wolfram Language code:

pndsol = ParametricNDSolve[{x'[t] == a x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-1) ^ t]}, x, {t, 0, 3}, {a}]

Wolfram Language code: Plot[(x[1] /. pndsol)[t], {t, 0, 3}]

NDSolveValue:

Wolfram Language code:

ndsol = NDSolveValue[{x'[t] == x[t] / 5 + 3, x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-3) ^ t]}, x, {t, 0, 3}]

Wolfram Language code: Plot[ndsol[t], {t, 0, 3}]

ParametricNDSolveValue:

Wolfram Language code:

pndsol = ParametricNDSolveValue[{x'[t] == x[t] / 5 + 3, x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-3) ^ t]}, x, {t, 0, 3}, {a}]

Wolfram Language code: Plot[pndsol[1][t], {t, 0, 3}]

DSolveValue:

Wolfram Language code: dsol = DSolveValue[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, x[t] -> x[t] + (-1) ^ t]}, x, {t, 0, 3}]

Wolfram Language code: Plot[dsol[t], {t, 0, 3}, Exclusions -> None]

推广和延伸 (2)

有条件的停止积分：

Wolfram Language code:

NDSolve[{x''[t] + .1 x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0, WhenEvent[x[t] == 0, Print[{t, x'[t]}];If[Abs[x'[t]] < .5, "StopIntegration"]]}, x, {t, ∞}]

有条件的去除一个事件：

Wolfram Language code:

Plot[Evaluate[{x[t], x'[t]} /. NDSolve[{x''[t] + .1 x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0, WhenEvent[x[t] == 0, Print[{t, x'[t]}];If[Abs[x'[t]] < .5, "RemoveEvent"]]}, x, {t, 0, 20}]], {t, 0, 20}]

选项 (11)

"DetectionMethod" (3)

默认事件检测的方法是 "Sign"，检验步骤间的符号变化：

Wolfram Language code: NDSolve[{x''[t] == -x[t], x[0] == 0, x'[0] == 1, WhenEvent[x[t] == 0, Print[t]]}, x, {t, 0, 10}]

上面等价于：

Wolfram Language code:

NDSolve[{x''[t] == -x[t], x[0] == 0, x'[0] == 1, WhenEvent[x[t] == 0, Print[t], "DetectionMethod" -> "Sign"]}, x, {t, 0, 10}];

Wolfram Language code: Plot[x[t] /. %, {t, 0, 10}]

如果步骤直接的交叉点数目是偶数，"Sign" 可能错过事件：

Wolfram Language code: NDSolve[{x'[t] == -1, x[0] == 1, WhenEvent[x[t] ^ 2 == .000001, Print[t], "DetectionMethod" -> "Sign"]}, x, {t, 0, 4}]

如果符号相同，"DerivativeSign" 将步长缩短为一半，但是导数符号改变：

Wolfram Language code:

NDSolve[{x'[t] == -1, x[0] == 1, WhenEvent[x[t] ^ 2 == .000001, Print[t], "DetectionMethod" -> "DerivativeSign"]}, x, {t, 0, 4}]

对于快速变化的事件函数，使用 "Interpolation" 方法定位事件：

Wolfram Language code:

{sol, events} = NDSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Sin[Pi x[t]] == 0, Sow[t], "DetectionMethod" -> "Interpolation"]}, x, {t, 0, 3}]//Reap;
{Length[First@events], Floor[Exp[3]]}

Wolfram Language code: Plot[{Sin[Pi x[t]]} /. sol, {t, 0, 3}]

"LocationMethod" (3)

默认定位方法 "Brent" 精确定位一个事件：

Wolfram Language code:

NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0, WhenEvent[y'[t] ≥ 0, {Print["Automatic", ": t = ", t, ", y'[t] = ", y'[t]], "StopIntegration"}]}, y, {t, 0, ∞}, Method -> "ExplicitRungeKutta"];

"LinearInterpolation" 方法会合理定位事件：

Wolfram Language code:

NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0, WhenEvent[y'[t] ≥ 0, {Print["LinearInterpolation", ": t = ", t, ", y'[t] = ", y'[t]], "StopIntegration"}, "LocationMethod" -> "LinearInterpolation"]}, y, {t, 0, ∞}, Method -> "ExplicitRungeKutta"];

如果高精度的事件定位不是必要的，则使用 "StepBegin" 或 "StepEnd"：

Wolfram Language code:

Map[
	NDSolve[{y''[t] + Sin[y[t]] == 0, y[0] == 3, y'[0] == 0, WhenEvent[y'[t] ≥ 0, {Print[#, ": t = ", t, ", y'[t] = ", y'[t]], "StopIntegration"}, "LocationMethod" -> #]}, y, {t, 0, ∞}, Method -> "ExplicitRungeKutta"]&, 
	{"StepBegin", "StepEnd"}
	];

"IntegrateEvent" (2)

当事件函数比解变化更快时，可能错过相交点：

Wolfram Language code:

{sol, {events, steps}} = Reap[NDSolveValue[{x'[t] == x[t], x[0] == 1, WhenEvent[Sin[Pi x[t]] == 0, Sow[t, "Events"]]}, x, {t, 0, 3}, StepMonitor :> Sow[t, "Steps"]]];

在与之间绘制具有时间步长和交叉点的事件函数：

Wolfram Language code:

Plot[Sin[π sol[t]], {t, 2, 3}, Epilog -> {{Red, Point[Thread[{events, 0}]]}, {Gray, Map[Line[{{#, -1}, {#, 1}}]&, steps]}}]

"IntegrateEvent"->True 时，步长会取得更小一些：

Wolfram Language code:

{sol, {events, steps}} = Reap[NDSolveValue[{x'[t] == x[t], x[0] == 1, WhenEvent[Sin[Pi x[t]] == 0, Sow[t, "Events"], "IntegrateEvent" -> True]}, x, {t, 0, 3}, StepMonitor :> Sow[t, "Steps"]]];
Plot[Sin[π sol[t]], {t, 2, 3}, Epilog -> {{Red, Point[Thread[{events, 0}]]}, {Gray, Map[Line[{{#, -1}, {#, 1}}]&, steps]}}]

定义一个函数，比较事件检测性能：

Wolfram Language code:

perf[T_][integrate_, detection_] := Flatten[AbsoluteTiming[Block[{e = 0}, NDSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Sin[Pi x[t]] == 0, e++, "IntegrateEvent" -> integrate, "DetectionMethod" -> detection]}, x, {t, 0, T}, MaxSteps -> ∞];{integrate, detection, Floor[Exp[T]] - e}]]]

"DetectionMethod"->"Interpolation" 时，默认情况下，事件进行鲁棒性积分：

Wolfram Language code: perf[10][Automatic, "Interpolation"]

比较性能：

Wolfram Language code:

TableForm[Flatten[Outer[perf[10], {True, False}, {"Interpolation", "DerivativeSign", "Sign"}], 1], TableHeadings -> {{}, {"Timing", "IntegrateEvent", "DetectionMethod", "Events missed"}}]

"Priority" (3)

当事件同时发生时，设置优先级别，该例是为每个整数：

Wolfram Language code:

NDSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1 / 2] == 0, Print["t = ", t, ", priority 2 event"], "Priority" -> 2], WhenEvent[Mod[t, 1] == 0, Print["t = ", t, ", priority 1 event"], "Priority" -> 1]}, x, {t, 0, 2}];

在任何有限优先级之后，但在优先级 Infinity 之前，将产生 Automatic 优先级：

Wolfram Language code:

NDSolve[{x'[t] == x[t], x[0] == 1, WhenEvent[Mod[t, 1] == 0, Print["t = ", t, ", Automatic priority event"]], WhenEvent[Mod[t, 1] == 0, Print["t = ", t, ", priority ", Infinity, " event"], "Priority" -> Infinity], 
	WhenEvent[Mod[t, 1] == 0, Print["t = ", t, ", priority 1 event"], "Priority" -> 1]}, {x}, {t, 0, 1}];

Sort 是按标准顺序对优先级排序：

Wolfram Language code:

plist = {Automatic, 2, 1, Infinity, 4};NDSolve[{x'[t] == 1, x[0] == 0, Map[WhenEvent[x[t] == 1, Print["Priority ", #], "Priority" -> #]&, plist]}, x, {t, 0, 2}];Sort[plist]

应用 (21)

反弹球 (3)

模拟一个球反弹下的步骤：

Wolfram Language code:

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

Wolfram Language code: Plot[Evaluate[{y[t], a[t]} /. sol], {t, 0, 8}, Filling -> {2 -> 0}]

绘制球的动能、势能和总能量：

Wolfram Language code:

kin[v_] := .5 v ^ 2;
pot[y_] := 9.8y;
energy[y_, v_] := kin[v] + pot[y];

Wolfram Language code: Plot[Evaluate[{kin[y'[t]], pot[y[t]], energy[y[t], y'[t]]} /. sol], {t, 0, 8}]

模拟20个反弹球同时从不同的高度下降：

Wolfram Language code:

With[{n = 20}, sol = NDSolve[Table[With[{i = i}, {y[i]''[t] == g[i][t], y[i][0] == i, y[i]'[0] == 0, g[i][0] == -9.81, WhenEvent[y[i][t] == 0, If[Abs[y[i]'[t]] > 10 ^ -6, y[i]'[t] -> -0.7 y[i]'[t], {y[i]'[t], g[i][t]} -> {0, 0}]]}], {i, Range[n]}], Array[y, n], {t, 0, 8}, DiscreteVariables -> Array[g, n]];
	Plot[Evaluate[Table[y[i][t], {i, Range[n]}] /. sol], {t, 0, 8}, PlotRange -> All]]

模拟一个球的正方形内的侧墙的影响下改变方向的运动：

Wolfram Language code:

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

Wolfram Language code: ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False]

具有无理初始速度，不再有周期解：

Wolfram Language code:

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

Wolfram Language code: ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False]

最优化 (1)

通过事件检测求解微分方程可用于查找函数在有限范围内的所有局部最小值和最大值. 考虑函数：

Wolfram Language code: f = Sin[10Pi x ^ 3] / (2x) + (x - 1 / 2) ^ 4;

设置搜索范围：

Wolfram Language code: xLimits = {-1.3, 1.3};

求解函数导数为 0 的微分方程并存储点：

Wolfram Language code:

{sol, points} = Reap@NDSolve[{y'[x] == D[f, x], y[xLimits[[1]]] == (f /. x -> xLimits[[1]]), WhenEvent[y'[x] == 0, Sow[{x, y[x]}]]}, y[x], {x, xLimits[[1]], xLimits[[2]]}];

绘制函数以及最小值和最大值：

Wolfram Language code:

Plot[f, {x, xLimits[[1]], xLimits[[2]]}, Epilog -> {PointSize[Medium], Red, Point[Flatten[points, 1]]}, PlotPoints -> 20]

交叉检测 (3)

标注参数曲线弧长度为1的点：

Wolfram Language code:

{sol, points} = Reap@
	NDSolve[{x'[t] == -y[t] - x[t] ^ 2, y'[t] == 2x[t] - y[t] ^ 3, p'[t] == Sqrt[Derivative[1][x][t]^2 + Derivative[1][y][t]^2], x[0] == y[0] == 1, p[0] == 0, WhenEvent[Mod[p[t], 1] == 0, Sow[t]]}, {x, y, p}, {t, 100}];

Wolfram Language code: ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 100}, Mesh -> points, MeshShading -> {Green, Orange}, PlotRange -> All]

标注参数解进入或离开环形区域的解：

Wolfram Language code:

wheninout = {WhenEvent[x[t] ^ 2  + y[t] ^ 2 == 10000, Sow[{x[t], y[t]}]], WhenEvent[x[t] ^ 2  + y[t] ^ 2 == 5000, Sow[{x[t], y[t]}]]};

Wolfram Language code:

{sol, points} = Reap@NDSolve[{x'[t] == -t y[t] + x[t], y'[t] == t x[t] + y[t], x[0]  == 1, y[0] == 0, wheninout}, {x[t], y[t]}, {t, 0, 10}];

Wolfram Language code:

Show[RegionPlot[5000 <= x ^ 2 + y ^ 2 <= 10000, {x, -110, 110}, {y, -110, 110}, Epilog -> {Red, PointSize[Medium], Point@@points}], ParametricPlot[{x[t], y[t]} /. sol, {t, 0, 5}]]

只标注参数解离开区域的点：

Wolfram Language code:

eqns = {x'[t] - y[t] == Cos[ϵ t], y'[t] + x[t] == 3Sin[ϵ t], x[0] == 0, y[0] == 0};
 
whenout = {WhenEvent[(x[t] / 3) ^ 2 + y[t] ^ 2 > 1, Sow[{x[t], y[t]}]]};
ϵ = .1;

Wolfram Language code: {sol, points} = Reap@NDSolve[Join[eqns, whenout], {x[t], y[t]}, {t, 0, 100}];

Wolfram Language code:

Show[ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100}, Epilog -> {PointSize[Medium], Red, Point@@points}], RegionPlot[(x / 3) ^ 2 + y ^ 2 ≤ 1, {x, -3, 3}, {y, -2, 2}, PlotStyle -> {Yellow, Opacity[.1]}]]

只标注参数解进入区域的点：

Wolfram Language code: whenin = {WhenEvent[(x[t] / 3) ^ 2 + y[t] ^ 2 < 1, Sow[{x[t], y[t]}]]};

Wolfram Language code: {sol, points} = Reap@NDSolve[Join[eqns, whenin], {x[t], y[t]}, {t, 0, 100}];

Wolfram Language code:

Show[ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100}, Epilog -> {PointSize[Medium], Red, Point@@points}], RegionPlot[(x / 3) ^ 2 + y ^ 2 ≤ 1, {x, -3, 3}, {y, -2, 2}, PlotStyle -> {Yellow, Opacity[.1]}]]

标注解进入或离开区域的点：

Wolfram Language code: wheninout = {WhenEvent[(x[t] / 3) ^ 2 + y[t] ^ 2 == 1, Sow[{x[t], y[t]}]]};

Wolfram Language code: {sol, points} = Reap@NDSolve[Join[eqns, whenin, whenout], {x[t], y[t]}, {t, 0, 100}];

Wolfram Language code:

Show[ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 100}, Epilog -> {PointSize[Medium], Red, Point@@points}], RegionPlot[(x / 3) ^ 2 + y ^ 2 ≤ 1, {x, -3, 3}, {y, -2, 2}, PlotStyle -> {Yellow, Opacity[.1]}]]

混合动力系统 (3)

每次线性振荡解穿过负轴，反应它穿过轴：

Wolfram Language code:

de[δ_] := {Subscript[x, 1]'[t] == Subscript[x, 2][t], Subscript[x, 2]'[t] == -Subscript[x, 1][t] + 2δ Subscript[x, 2][t] + 1};
ic = {Subscript[x, 1][0] == .5, Subscript[x, 2][0] == .5};

Wolfram Language code:

{sol, {points}} = NDSolve[{de[0.1], ic, WhenEvent[And[Subscript[x, 2][t] == 0, Subscript[x, 1][t] < 0], {Subscript[x, 1][t] -> -Subscript[x, 1][t], Sow[-Subscript[x, 1][t]]}]}, {Subscript[x, 1], Subscript[x, 2]}, {t, 0, 1000}]//Reap;

该复位振荡器的解展现出混沌的行为：

Wolfram Language code:

Map[ParametricPlot[{Subscript[x, 1][t], Subscript[x, 2][t]} /. sol, {t, 0, #}, PlotPoints -> 500, PlotRange -> {{-1, 3}, {-2, 2}}]&, {10, 100, 400}]

在具有反射点的直方图的负轴上绘制解：

Wolfram Language code:

Grid[{{Histogram[points, {-1, 0, .05}]}, {ParametricPlot[{Subscript[x, 1][t], Subscript[x, 2][t]} /. sol, {t, 0, 1000}, PlotRange -> {{-1, 0}, {-0.3, 0}}, PlotPoints -> 2000, Axes -> None]}}]

模拟具有正弦压迫的自由度为1的效应振荡器（impact oscillator）：

Wolfram Language code:

sol[ω_] := With[{r = 0.95}, 
	NDSolve[{u'[t] == v[t], v'[t] == -u[t] + Cos[ω t], u[0] == 0, v[0] == 0, WhenEvent[u[t] == 0, v[t] -> -r v[t]]}, {u, v}, {t, 0, 300}]];

Wolfram Language code:

Table[ParametricPlot[Evaluate[{u[t], v[t]} /. sol[ω]], {t, 200, 300}, PlotPoints -> 200, PlotRange -> {{-.1, 1}, {-1, 1}}], {ω, {3, 2.76, 2.9}}]

模拟一个在固定时间间隔被敲打的阻尼振荡器：

Wolfram Language code:

system = {x''[t] + .1x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0};
control = WhenEvent[Mod[t, 1], x'[t] -> x'[t] + 1];

Wolfram Language code: sol = NDSolve[{system, control}, x, {t, 0, 100}]

Wolfram Language code: Plot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 100}, PlotPoints -> 120]

轨迹最终安定为一个连贯的轨道：

Wolfram Language code: ParametricPlot[{x[t], x'[t]} /. sol, {t, 90, 100}, AspectRatio -> 1]

使用同样的振荡器，但是在固定的时间间隔进行随机的敲打：

Wolfram Language code:

system = {x''[t] + .1x'[t] + x[t] == 0, x[0] == 1, x'[0] == 0};
control = WhenEvent[Mod[t, 1], x'[t] -> x'[t] + 0.05RandomVariate[NormalDistribution[]]];

Wolfram Language code: sol = NDSolve[{system, control}, x, {t, 0, 100}]

Wolfram Language code: Plot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 100}, PlotPoints -> 120]

Wolfram Language code: ParametricPlot[{x[t], x'[t]} /. sol, {t, 0, 100}, AspectRatio -> 1]

摩擦模型 (1)

对使用弹簧固定在墙壁上的方块，使用摩擦力的不同模型，包括粘性、库伦力、Stribeck 和静态. 比较不同模型的位置和速度：

Wolfram Language code:

m = l = 1.;
beltv[t_] = .1;
spring[x_] = 1000.(l - x);

方块的牛顿方程：

Wolfram Language code: sys := {m x''[t] == spring[x[t]] + friction[x'[t]], x[0] == 1, x'[0] == 0};

粘性摩擦力 与相对速度成正比：

Wolfram Language code:

viscous[v_] := -30.(v - beltv[t]);
friction[v_] := viscous[v];

Wolfram Language code: {pos1, vel1} = NDSolveValue[sys, {x[t], x'[t]}, {t, 0, 2}];

方块在弹簧的自然长度为 1 之上稳定：

Wolfram Language code: {Plot[{pos1}, {t, 0, 1}, PlotRange -> All], Plot[{vel1, beltv[t]}, {t, 0, 1}, PlotRange -> All]}

库伦摩擦力 与相对速度的符号成正比：

Wolfram Language code:

coulomb[v_] := -25.Sign[v - beltv[t]];
friction[v_] := viscous[v] + coulomb[v];

Wolfram Language code: {pos2, vel2} = NDSolveValue[sys, {x[t], x'[t]}, {t, 0, 20}];

方块随着传送带移动，直至弹簧力足够大：

Wolfram Language code: {Plot[{pos2}, {t, 0, 1}, PlotRange -> All], Plot[{vel2, beltv[t]}, {t, 0, 1}, PlotRange -> All]}

Stribeck 摩擦力 是一种修正后的库伦摩擦力 $F_(str)=gamma sgn(v) e^(-2 TemplateBox[{v}, Abs])$ :

Wolfram Language code:

stribeck[v_] := -.3Sign[v]Exp[-2Abs[v]];
friction[v_] := viscous[v] + coulomb[v] + stribeck[v];

Wolfram Language code: {pos3, vel3} = NDSolveValue[sys, {x[t], x'[t]}, {t, 0, 20}];

低速上的变化略有减少：

Wolfram Language code: {Plot[{pos3}, {t, 0, 1}, PlotRange -> All], Plot[{vel3, beltv[t]}, {t, 0, 1}, PlotRange -> All]}

静摩擦力 将方块保持在原地直至取决于表面粗糙度，弹簧力超过一定值 μ. 当方块被卡住时，使用设为 1 的离散变量 stuck，否则使用 0：

Wolfram Language code: friction[v_] := viscous[v] + coulomb[v] + stribeck[v];

Wolfram Language code: staticsys := {x''[t] == If[stuck[t] == 1, beltv'[t], spring[x[t]] + friction[x'[t]]], x[0] == 1, x'[0] == 0};

检查弹簧力是否小于 μ，并且该方块是否相对于滑动带没有移动：

Wolfram Language code:

μ = 100;
stick = WhenEvent[x'[t] == beltv[t], stuck[t] -> Boole[spring[x[t]]^2 < μ^2]];
slip = WhenEvent[spring[x[t]]^2 > μ^2, stuck[t] -> 0];

Wolfram Language code:

{pos4, vel4} = NDSolveValue[{staticsys, stick, slip, stuck[0] == 0}, {x[t], x'[t]}, {t, 0, 8}, DiscreteVariables -> stuck[t]];

方块持续黏在传送带上，然后由于弹簧力滑走：

Wolfram Language code: {Plot[{pos4}, {t, 0, 2}, PlotRange -> All], Plot[{vel4, beltv[t]}, {t, 0, 2}, PlotPoints -> 100, PlotRange -> All]}

比较不同模型：

Wolfram Language code:

Plot[{pos1, pos2, pos3, pos4}, {t, 0, 2}, PlotRange -> {0.9, 1.1}, PlotLegends -> {"viscous", "viscous+Coulomb", "viscous+Coulomb+Stribeck", "viscous+Coulomb+Stribeck+static"}]

电力电子学 (4)

对由四个二极管和一个电容构成的 AC-to-DC 全波整流器建模，用于平滑输出：

输入电压 vi 和整流后的电压 vr：

Wolfram Language code: vi[t_] = Sin[t] ;

Wolfram Language code: vr[t_] = Abs[vi[t]]//PiecewiseExpand[#, t∈Reals]&;

Wolfram Language code: Plot[{vi[t], vr[t]}, {t, 0, 16}]

当 vr 增加时，vo[t]=vr[t]；电容器充电，并将电流供给到负载. 当 vr 开始降低 (vr'[t]<0)，电容器通过负载放电，输出电压服从 vo'[t]=vo[t]/c r ：

Wolfram Language code: system = {vo'[t] == charging[t]vr'[t] - (1 - charging[t])vo[t] / (r c), vo[0] == 0};

Wolfram Language code:

control = {charging[0] == 1, WhenEvent[vr'[t] < 0, charging[t] -> 0], WhenEvent[And[vo[t] < vr[t], charging[t] == 0, vr'[t] > 0], charging[t] -> 1]};

对系统进行仿真：

Wolfram Language code: pars = {r -> 0.1, c -> 1000};

Wolfram Language code: sol = NDSolve[{system, control} /. pars, vo, {t, 0, 16}, DiscreteVariables -> charging∈{0, 1}];

绘制电容平滑处理后的输出电压：

Wolfram Language code: Plot[{vo[t] /. sol, vr[t]}, {t, 0, 16}]

对使用脉冲宽度调制反馈控制 q[t] 的 DC-to-DC 降压转换器建模，该转换器将输入电压 vi 转换到所需的输出电压 vd：

使用基尔霍夫定理获得上述电路的模型:

Wolfram Language code: system = {vo'[t] == -vo[t] / (r c) + i[t] / c , i'[t] == -vo[t] / l + q[t]vi / l, vo[0] == 0, i[0] == 0};

控制信号 q[t] 将在占每个周期 τ 的 vd/vi 时间段内接通晶体管：

Wolfram Language code: control = {q[0] == 1, WhenEvent[Mod[t, τ] == (vd / vi)τ, q[t] -> 0], WhenEvent[Mod[t, τ] == 0, q[t] -> 1]};

从一个较高的电压 vi=24 降到较低的电压 vd=16:

Wolfram Language code: pars = {vi -> 24, vd -> 16, r -> 22, l -> 2 10 ^ -2, c -> 4 10 ^ -6, τ -> 5 10 ^ -5};

Wolfram Language code: sol = NDSolve[{system, control} /. pars, {vo, i, q}, {t, 0, .01}, DiscreteVariables -> q, MaxSteps -> Infinity];

Wolfram Language code: {Plot[{vo[t] /. sol, vd /. pars}, {t, 0, .01}, PlotRange -> All], Plot[q[t] /. sol, {t, 0, .0004}]}

对使用脉冲宽度调制反馈控制 q[t] 的 DC-to-DC 升压转换器建模，该转换器将输入电压 vi 转换到所需的输出电压 vd：

使用基尔霍夫定理获得上述电路的模型：

Wolfram Language code: system = {vo'[t] == q[t]i[t] / c - vo[t] / (r c), i'[t] == -q[t]vo[t] / l + vi / l, vo[0] == 0, i[0] == 0};

控制信号 q[t] 将在占每个周期 τ 的 vi/vd 时间段内接通晶体管：

Wolfram Language code: control = {q[0] == 1, WhenEvent[Mod[t, τ] == (vi / vd)τ, q[t] -> 0], WhenEvent[Mod[t, τ] == 0, q[t] -> 1]};

从较低的电压 vi=24 升到较高的电压 vd=36:

Wolfram Language code: pars = {vi -> 24, vd -> 36, r -> 22, l -> 2 10 ^ -1, c -> 1 10 ^ -4, τ -> 5 10 ^ -5};

Wolfram Language code: sol = NDSolve[{system, control} /. pars, {vo, i, q}, {t, 0, .2}, DiscreteVariables -> q, MaxSteps -> Infinity];

Wolfram Language code: {Plot[{vo[t] /. sol, vd /. pars}, {t, 0, .2}, PlotRange -> All], Plot[q[t] /. sol, {t, 0, .0004}]}

对从输入电压电平 vi 到所需输出电压电平 vo 的 DC-to-DC 升压转换器 使用脉冲宽度调制反馈控制 q[t] 建模：

使用基尔霍夫定理获得上述电路的模型：

Wolfram Language code: system = {i'[t] == q[t]vi / l - (1 - q[t])vo[t] / l , vo'[t] == i[t] / c - vo[t] / (r c), vo[0] == 0, i[0] == 0};

控制信号 q[t] 将晶体管开启每个周期的 vd/(vi+vd) 时间：

Wolfram Language code: control = {q[0] == 1, WhenEvent[Mod[t, τ] == (vd / (vd + vi))τ, q[t] -> 0], WhenEvent[Mod[t, τ] == 0, q[t] -> 1]};

从较低的电压 vi=24 升到较高的电压 vd=36:

Wolfram Language code: pars = {vi -> 24, vd -> 36, r -> 22, l -> 2 10 ^ -2, c -> 4 10 ^ -6, τ -> 5 10 ^ -5};

Wolfram Language code: sol = NDSolve[{system, control} /. pars, {vo, i, q}, {t, 0, .2}, DiscreteVariables -> q, MaxSteps -> Infinity];

Wolfram Language code: {Plot[{vo[t] /. sol, vd /. pars}, {t, 0, .02}, PlotRange -> All], Plot[q[t] /. sol, {t, 0, .0004}]}

从较高的电压 vi=24 降到较低的电压 vd=16:

Wolfram Language code: pars = {vi -> 24, vd -> 16, r -> 22, l -> 2 10 ^ -2, c -> 4 10 ^ -6, τ -> 5 10 ^ -5};

Wolfram Language code: sol = NDSolve[{system, control} /. pars, {vo, i, q}, {t, 0, .2}, DiscreteVariables -> q, MaxSteps -> Infinity];

Wolfram Language code: {Plot[{vo[t] /. sol, vd /. pars}, {t, 0, .02}, PlotRange -> All], Plot[q[t] /. sol, {t, 0, .0004}]}

数字控制 (5)

仿真用离散时间控制器稳定的系统：

Wolfram Language code:

system = {y'[t] == y[t] + u[t], y[0] == 1, u[0] == 0};
control = WhenEvent[Mod[t, 1.0], u[t] -> -2y[t]];

仿真和可视化：

Wolfram Language code: NDSolve[{system, control}, {y, u}, {t, 0, 20}, DiscreteVariables -> u];

Wolfram Language code: Plot[Evaluate[{u[t], y[t]} /. %], {t, 0, 10}, PlotRange -> {-6, 6}]

使用无振 (dead-beat) 离散时间控制器控制双重积分器：

Wolfram Language code: system = {x''[t] == u[t], x[0] == 1, x'[0] == 1, u[0] == -7};

使用无振（dead-beat）数字反馈控制器：

Wolfram Language code:

τ = 0.5;
control = WhenEvent[Mod[t, τ], u[t] -> -(1 / τ ^ 2)x[t] - (3 / (2τ))x'[t]];

仿真和可视化：

Wolfram Language code: NDSolve[{system, control}, {x, x', u}, {t, 0, 2}, DiscreteVariables -> u];

Wolfram Language code: Plot[Evaluate[{x[t], x'[t], u[t]} /. %], {t, 0, 2}, PlotRange -> All]

用 1 公斤质体模拟移动物体的位置：

Wolfram Language code: system = {x''[t] == u[t], x'[0] == x[0] == 0, u[0] == 1};

用采样比例-微分（PD）控制器保持位置常量：

Wolfram Language code:

kp = 1; td = 1;τ = 1; xref = 1;
control = WhenEvent[Mod[t, τ], u[t] -> kp(xref - x[t] - td x'[t])];

仿真和可视化：

Wolfram Language code: NDSolve[{system, control}, {x, u}, {t, 0, 12}, DiscreteVariables -> u];

Wolfram Language code: Plot[Evaluate[{xref, x[t], u[t]} /. %], {t, 0, 12}]

设计具有状态方程，和输出的机械设备（plant）的离散时间控制器，仿真闭环系统：

Wolfram Language code: plant = {Subscript[x, 1]'[t] == 2 Subscript[x, 1][t] + u[t], Subscript[x, 2]'[t] == 3Subscript[x, 2][t] + u[t]};

Wolfram Language code:

pss = StateSpaceModel[plant, {Subscript[x, 1][t], Subscript[x, 2][t]}, {u[t]}, Subscript[x, 1][t] + Subscript[x, 2][t], t]

通过使用状态估计器和状态反馈设计连续时间的控制器：

Wolfram Language code: ctreg = EstimatorRegulator[pss, {EstimatorGains[pss, {-12, -14}], StateFeedbackGains[pss, {-6, -7}]}]

用采样时间采样获取一个离散时间控制器：

Wolfram Language code: τ = 0.0001;

Wolfram Language code: dtreg = ToDiscreteTimeModel[ctreg, τ, Method -> "ZeroOrderHold"]

使用离散时间控制器:

Wolfram Language code: {a, b, c, d} = dtreg//Normal;

Wolfram Language code:

controller = Sequence[WhenEvent[Mod[t, τ], {Subscript[z, 1][t], Subscript[z, 2][t], Overscript[u,  ^ ][t]} -> 
	Evaluate@Join[a.{Subscript[z, 1][t], Subscript[z, 2][t]} + b.{u[t], Subscript[x, 1][t] + Subscript[x, 2][t]}, c.{Subscript[z, 1][t], Subscript[z, 2][t]} + d.{u[t], Subscript[x, 1][t] + Subscript[x, 2][t]}]], u[t] == -Overscript[u,  ^ ][t]];

仿真具有控制器和 plant 的闭环系统：

Wolfram Language code: init = {Subscript[x, 1][0] == 0.1, Subscript[x, 2][0] == 0.1, Subscript[z, 1][0] == 0., Subscript[z, 2][0] == 0};

Wolfram Language code:

sol = NDSolve[{plant, controller, init}, {Subscript[x, 1], Subscript[x, 2], u}, {t, 0, 2}, DiscreteVariables -> {Overscript[u,  ^ ], Subscript[z, 1], Subscript[z, 2]}];

绘制结果：

Wolfram Language code: Plot[Evaluate[{Subscript[x, 1][t] + Subscript[x, 2][t], u[t]} /. sol], {t, 0, 2}, PlotRange -> All]

模拟一个固定在移动车上位置为 x，角度为 θ，承受水平力为 f 的钟摆：

Wolfram Language code:

system = {2 f[t] + Cos[θ[t]] Derivative[1][θ][t]^2 + Sin[θ[t]] Derivative[2][θ][t] == 4 Derivative[2][x][t], 2 Cos[θ[t]] - 2 Sin[θ[t]] Derivative[2][x][t] + Derivative[2][θ][t] == 0};

设计一个稳定关于不稳定垂直固定点钟摆的控制器：

Wolfram Language code: ssm = StateSpaceModel[system, {{x[t], 0}, {x'[t], 0}, {θ[t], (π/2)}, {θ'[t], 0}}, f[t], {}, t]

Wolfram Language code: gains = DiscreteLQRegulatorGains[N[ssm], {DiagonalMatrix[{1, 10, 10, 100}], {{1}}}, .1]//First

利用 WhenEvent 在固定采样时间应用反馈：

Wolfram Language code: control = {WhenEvent[Mod[t, .1], df[t] -> -gains.{x[t], x'[t], θ[t] - (π/2), θ'[t]}]};

仿真该系统：

Wolfram Language code:

sol = NDSolve[Join[system /. f[t] -> df[t], control, {df[0] == 0, x[0] == x'[0] == θ'[0] == 0, θ[0] == (π/2) - 1.}], {x, θ, df}, {t, 0, 15}, DiscreteVariables -> df[t]];

在倒置位置（）钟摆的角度快速稳定：

Wolfram Language code: Plot[Evaluate[{θ[t], df[t]} /. sol], {t, 0, 8}, PlotRange -> All]

PDE 模型 (1)

模拟温控器控制下一个房间中的热生成，该房间有三面绝热墙和一面受外界温度影响的玻璃墙：

Wolfram Language code:

Ω = Rectangle[{0, 0}, {3 / 2, 1}];
outsideTemp[t_] := 15 + 10 * Sin[2π t / 24];
kd = 0.78;Ld = 0.05;
Γ = NeumannValue[Ld / kd * (outsideTemp[t] - u[t, x, y]), {x == 0}];

在事件发生时，加热器的负载会增大或减小. 系数 20 指定过渡的斜率：

Wolfram Language code:

heater[upQ_, t_, tEvent_] := Evaluate[With[{heaterLoad = 26, slope = 20}, If[upQ == 1, Min[20 * Max[(t - tEvent), 0], 1], 1 - Min[20 * Max[(t - tEvent - 1 / 8), 0], 1]] * heaterLoad]]

可视化加热器在任意事件时间 10 处的情况：

Wolfram Language code:

With[{tEvent = 10}, 
	Plot[{heater[1, t, tEvent], heater[0, t, tEvent]}, {t, 9, 11}, PlotLegends -> {"up", "down"}]]

注意，up 和 down 阶段有重叠. 这可以防止它们在两种状态之间连续跳转，就像真正的控制器一样.

PDE 模拟空气中的热扩散，在圆内生成热，通过玻璃窗散热：

Wolfram Language code:

ρ = 1.225;Cp = 1005.4;
With[{heating = heater[a[t], t, eventT[t]]}, 
	pde = D[u[t, x, y], t]  - ρ * Cp * Laplacian[u[t, x, y], {x, y}] == If[(x - 1 / 2) ^ 2 + (y - 1 / 2) ^ 2 ≤ (2 / 10) ^ 2, heating, 0] + Γ];

如果位于处的温控器测量到的温度高于/低于触发温度，且离散变量发生变化，则关闭/打开加热器：

Wolfram Language code:

triggerLow = 18;triggerHigh = 20;
events = {a[0] == 1, eventT[0] == 0, WhenEvent[u[t, 1.25, .25] < triggerLow, {eventT[t], a[t]} -> {If[a[t] == 0, t, eventT[t]], 1}], WhenEvent[u[t, 1.25, .25] > triggerHigh, {eventT[t], a[t]} -> {If[a[t] == 1, t, eventT[t]], 0}]};

初始状态等于外部温度，查看 PDE 的时间积分：

Wolfram Language code:

eqn = {pde, u[0, x, y] == outsideTemp[0], events};
res = Monitor[NDSolveValue[eqn, {u, a}, {t, 0, 2 * 24}, {x, y}∈Ω, DiscreteVariables -> {eventT[t], a[t]}, EvaluationMonitor :> (monitor = Row[{"t = ", CForm[t]}])], monitor]

可视化温控器测量的温度、外部温度和加热器的触发. 加热器处于开启状态时显示蓝色背景：

Wolfram Language code:

hp = Plot[ 25res[[2]][t], {t, 0, 2 * 24}, Filling -> Bottom, PlotStyle -> None];Show[Plot[{res[[1]][t, 1.25, .25], outsideTemp[t], 18, 20}, {t, 0, 2 * 24}], hp]

查看积分过程中的时间步长，注意 NDSolveValue 如何在事件期间根据需要调整时间步长：

Wolfram Language code: ListPlot[res[[1]]["Coordinates"][[1]], AxesLabel -> {"Steps", "Time"}]

属性和关系 (2)

NDSolve 会尝试自动设置 WhenEvent 事件来处理中的不连续点：

Wolfram Language code: NDSolve[{x'[t] == Sign[1 - x[t]], x[0] == 2}, x, {t, 0, 2}];

Wolfram Language code: Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}]

在自动不连续处理后，下面等价问题将被解决：

Wolfram Language code:

NDSolve[{Derivative[1][x][t] == v[t], WhenEvent[1 == x[t], v[t] -> "DiscontinuitySignature"], x[0] == 2, v[0] == 1}, x, {t, 0, 2}, DiscreteVariables -> {Element[v, {-1, 0, 1}]}];

如果是一个黑盒（black-box）函数，NDSolve 不能自动处理不连续性：

Wolfram Language code: blackbox[x_ ? NumberQ] := Sign[1 - x];

Wolfram Language code: NDSolve[{x'[t] == blackbox[x[t]], x[0] == 2}, x, {t, 0, 2}];

这种情况，可以通过添加 WhenEvent 正确处理不连续性：

Wolfram Language code:

NDSolve[{x'[t] == v[t], x[0] == 2, v[0] == 1, WhenEvent[1 == x[t], v[t] -> "DiscontinuitySignature"]}, {x, v}, {t, 0, 2}, DiscreteVariables -> {Element[v, {-1, 0, 1}]}];

Wolfram Language code: Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}]

或者，使用 "CrossSlidingDiscontinuity" 穿越不连续性：

Wolfram Language code: NDSolve[{x'[t] == blackbox[x[t]], x[0] == 2, WhenEvent[1 - x[t], "CrossSlidingDiscontinuity"]}, x, {t, 0, 2}];

Wolfram Language code: Plot[{x[t] /. %, x'[t] /. %}, {t, 0, 2}]

可能存在的问题 (6)

事件诠释 (1)

显式给出条件的事件与单个谓词略有不同的诠释：

Wolfram Language code:

pred[t_ ? NumberQ, x_ ? NumberQ] := t > .5 && x > 0;
xsol = NDSolveValue[{x'[t] == 2 + x[t], x[0] == -1, WhenEvent[pred[t, x[t]], Print[{t, x[t]}]]}, x, {t, 0, 1}];

WhenEvent 不能找到显式根函数，因此当 pred 变成 True 时事件发生：

Wolfram Language code: Plot[{xsol[t], Boole[pred[t, xsol[t]]]}, {t, 0, 1}]

当给定显式函数 ()，事件被认为根，但条件是谓词 () 在时为 False，因此从不会被触发：

Wolfram Language code: xsol = NDSolveValue[{x'[t] == 2 + x[t], x[0] == -1, WhenEvent[t > .5 && x[t] > 0, Print[{t, x[t]}]]}, x, {t, 0, 1}];

显示解的踪迹和事件位置，背景为绿，当条件为 True，当条件为 False 时，背景为红色：

Wolfram Language code: Show[Plot[xsol[t], {t, 0, 1}], Graphics[...]]

当给定显式函数 ()，事件被作为根：

Wolfram Language code: xsol = NDSolveValue[{x'[t] == 2 + x[t], x[0] == -1, WhenEvent[x[t] > 0 && t > .5, Print[{t, x[t]}]]}, x, {t, 0, 1}];

显示解的踪迹和事件位置，背景为绿，当条件为 True，当条件为 False 时，背景为红色：

Wolfram Language code: Show[Plot[xsol[t], {t, 0, 1}], Graphics[...]]

事件检测 (2)

使用默认的 "DetectionMethod"，可能会错过一些事件：

Wolfram Language code: NDSolve[{y'[t] == -1, y[0] == 1, WhenEvent[y[t] ^ 2 == .000001, Print[t]]}, y, {t, 0, 4}]

尝试 "DerivativeSign" 或 "Interpolation" 获取更好的事件检测：

Wolfram Language code:

NDSolve[{y'[t] == -1, y[0] == 1, WhenEvent[y[t] ^ 2 == .000001, Print[t], "DetectionMethod" -> "DerivativeSign"]}, y, {t, 0, 4}]

Wolfram Language code:

NDSolve[{y'[t] == -1, y[0] == 1, WhenEvent[y[t] ^ 2 == .000001, Print[t], "DetectionMethod" -> "Interpolation"]}, y, {t, 0, 4}]

对应于事件顺序的根的事件可能会错过：

Wolfram Language code: NDSolve[{y'[t] == y[t], y[0] == 1, WhenEvent[(y[t] - 2) ^ 2 == 0, Print[t]]}, y, {t, 0, 1}]

相反，改写奇数根：

Wolfram Language code: NDSolve[{y'[t] == y[t], y[0] == 1, WhenEvent[(y[t] - 2) == 0, Print[t]]}, y, {t, 0, 1}]

任意接近事件 (1)

有些系统允许在有限的时间跨度中有无限数量的事件：

Wolfram Language code:

eqns = {y''[t] == -10, y[0] == 1, y'[0] == 0};
whenevent = {WhenEvent[y[t] == 0, {Sow[t], y'[t] -> -(7 / 10)y'[t]}]};
{sol, {points}} = NDSolve[Join[eqns, whenevent], y, {t, 0, 4}]//Reap;

最终，一个事件没有被检测到：

Wolfram Language code: Plot[y[t] /. sol, {t, 0, 2.7}, Epilog -> {Red, Point[({#, 0}&) /@ points]}]

Wolfram Language code: ListPlot[RealExponent[Differences[points]]]

为定位事件收紧公差，改善检测：

Wolfram Language code:

whenevent = {WhenEvent[y[t] == 0, {Sow[t], y'[t] -> -(7 / 10)y'[t]}, "LocationMethod" -> {"Brent", AccuracyGoal -> 12, PrecisionGoal -> 12}]};

Wolfram Language code: {sol, {times}} = NDSolve[Join[eqns, whenevent], y, {t, 0, 4}]//Reap;

Wolfram Language code: ListPlot[RealExponent[Differences[times]]]

如果比这更快的事件将不会被检测到：

Wolfram Language code: zp = Last[times];mindelta = Block[{delta = zp}, While[zp + delta > zp, delta = delta / 2]; delta]

通过插值检测，事件检测会进一步改善：

Wolfram Language code: whenevent = {WhenEvent[y[t] == 0, {Sow[t], y'[t] -> -(7 / 10)y'[t]}, "DetectionMethod" -> "Interpolation"]};

Wolfram Language code: {sol, {times}} = NDSolve[Join[eqns, whenevent], y, {t, 0, 4}]//Reap;

最后两个事件的间隔很接近这个极限：

Wolfram Language code: Last[Differences[times]]

在事件错过前自动停止积分：

Wolfram Language code:

mindelta = Block[{delta = 4.}, While[zp + delta > zp, delta = delta / 2]; delta];
	wheneventstop = {WhenEvent[y[t] == 0, {If[t - tlast ≤ 2mindelta, tend = t;"StopIntegration", tlast = Sow[t]], y'[t] -> -(7 / 10)y'[t]}, "DetectionMethod" -> "Interpolation"]};

Wolfram Language code:

Block[{tlast = -1}, {sol, {times}} = NDSolve[{eqns, wheneventstop}, y, {t, 0, 4}]//Reap;]
ListPlot[RealExponent[Differences[points]]]

事件在积分停止点被检测到：

Wolfram Language code: Plot[y[t] /. sol, {t, tend - 10^-12, tend}]

增加工作精度允许在更精细的时间尺度中检测事件：

Wolfram Language code: {sol, {times}} = NDSolve[{eqns, whenevent}, y, {t, 0, 10}, WorkingPrecision -> 50]//Reap;

Wolfram Language code: ListPlot[RealExponent[Differences[times]]]

事件行为 (2)

在 ODE 中不能重置最高阶导数：

Wolfram Language code: xsol = NDSolveValue[{x'[t] == -3x[t], x[0] == 1 / 3, WhenEvent[x[t] == 1 / 9, x'[t] -> -1]}, x, {t, 0, 2}]

Wolfram Language code: Plot[{xsol[t], xsol'[t]}, {t, 0, 2}]

但作为 DAE 求解，可以重置最高阶导数：

Wolfram Language code:

xsol = NDSolveValue[{x'[t] == -3x[t], x[0] == 1 / 3, WhenEvent[x[t] == 1 / 9, x'[t] -> -1]}, x, {t, 0, 2}, Method -> {"EquationSimplification" -> "Residual"}]

通过重新初始化，求解对应的值，因此的残差为0：

Wolfram Language code: Plot[{xsol[t], xsol'[t]}, {t, 0, 2}]

对于有序事件行为，变量是按序被修改的：

Wolfram Language code:

sol = NDSolve[{y'[t] == 0, x'[t] == 0, x[0] == 1, y[0] == 2, WhenEvent[t == .5, {x[t] -> y[t], y[t] -> x[t]}]}, {x, y}, {t, 0, 1}];

Wolfram Language code: Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 1}, PlotRange -> All, PlotStyle -> {Blue, Red}, AxesOrigin -> {0, 0}]

若要交换变量值，使用同步事件：

Wolfram Language code:

sol = NDSolve[{y'[t] == 0, x'[t] == 0, x[0] == 1, y[0] == 2, WhenEvent[t == .5, {x[t], y[t]} -> {y[t], x[t]}]}, {x, y}, {t, 0, 1}];

Wolfram Language code: Plot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 1}, PlotRange -> All, PlotStyle -> {Blue, Red}, AxesOrigin -> {0, 0}]

互动范例 (2)

用有理和无理初始值比较一个在盒中的球：

Wolfram Language code:

ball2d[r_ ? NumericQ] := NDSolveValue[{
	x'[t] == vx[t], WhenEvent[x[t] ^ 2 == 1, vx[t] -> -vx[t]], 
	y'[t] == vy[t], WhenEvent[y[t] ^ 2 == 1, vy[t] -> -vy[t]], 
	x[0] == 0, y[0] == 0, vx[0] == 1, vy[0] == r
	}, {x[t], y[t]}, {t, 0, 200}, DiscreteVariables -> {vx, vy}];

Wolfram Language code: {sol1, sol2} = {ball2d[4 / 3], ball2d[Sqrt[2]]};

Wolfram Language code:

Manipulate[
	GraphicsRow[ParametricPlot[#, {t, 0, τ}, Frame -> True, FrameTicks -> None, PlotRange -> 1, Axes -> False, PlotPoints -> 500, PlotStyle -> Thickness[0.02], Epilog -> {PointSize[.05], Point[# /. (t -> τ)]}]& /@ {sol1, sol2}], {{τ, 11.9}, 0.01, 198}, SaveDefinitions -> True]

用有理和无理初始值比较一个在盒中的球：

Wolfram Language code:

ball3d[r1_ ? NumericQ, r2_ ? NumericQ] := NDSolveValue[{x'[t] == a[t], y'[t] == b[t], z'[t] == c[t], x[0] == 0, y[0] == 0, z[0] == 0, a[0] == 1, b[0] == r1, c[0] == r2, WhenEvent[x[t] ^ 2 - 1, a[t] -> -a[t]], WhenEvent[y[t] ^ 2 - 1, b[t] -> -b[t]], WhenEvent[z[t] ^ 2 - 1, c[t] -> -c[t]]}, {x[t], y[t], z[t]}, {t, 0, 200}, DiscreteVariables -> {a, b, c}];

Wolfram Language code: {sol1, sol2} = {ball3d[1 / 3, 2 / 3], ball3d[Sqrt[2], Sqrt[3]]};

Wolfram Language code:

Manipulate[
	GraphicsRow[Show[ParametricPlot3D[#, {t, 0, τ}, Boxed -> True, PlotRange -> 1, Axes -> False, PlotPoints -> 500], Graphics3D[{Blue, Sphere[# /. (t -> τ), 1 / 8] /. (t -> τ)}]]& /@ {sol1, sol2}], {{τ, 11.5}, 0.01, 200, 1}, SaveDefinitions -> True]

Top

	f0	解的变量函数 f 跨过零
	f>0	从下面 f 跨过零
	f<0	从上面 f 跨过零
	Mod[t,]0	以固定间隔，时间变量采样

	"StopIntegration"	在停止积分，返回解
	"RestartIntegration"	在重新开始积分
	"CrossDiscontinuity"	积分到，推算并在重新开始
	"CrossSlidingDiscontinuity"	积分到，推算并检查滑动模式条件，在重新开始
	"RemoveEvent"	去除事件
	y[t]->val	把状态变量 y 改为 val
	d[t]->"DiscontinuitySignature"	改变不连续签名变量 d

	Automatic	自动决定检测方法
	"Sign"	使用一个符号变换
	"DerivativeSign"	使用符号变换和时间导数
	"Interpolation"	使用稠密解输出的插值

	Automatic	自动确定定位方法
	"StepBegin"	在每个步骤开始前定位事件
	"StepEnd"	在每个步骤结束时定位事件
	"LinearInterpolation"	在步骤中线性插值
	"Brent"	在一个步骤中使用 Brent 的根定位

	Automatic	自动决定
	False	不积分事件函数
	True	积分事件函数

WhenEvent

更多信息和选项

范例

基本范例 (1)

范围 (16)

事件 (6)

行为 (9)

求解器 (1)

推广和延伸 (2)

选项 (11)

"DetectionMethod" (3)

"LocationMethod" (3)

"Priority" (3)

应用 (21)

反弹球 (3)

最优化 (1)

交叉检测 (3)

混合动力系统 (3)

摩擦模型 (1)

电力电子学 (4)

数字控制 (5)

PDE 模型 (1)

属性和关系 (2)

可能存在的问题 (6)

事件诠释 (1)

事件检测 (2)

任意接近事件 (1)

事件行为 (2)

互动范例 (2)

参见

相关指南

历史

文本

CMS

APA

BibTeX

BibLaTeX