请问各位大佬，mathematica在数值求解微分方程
mathematica吧
全部回复
仅看楼主
level 1
小张同学hello3 楼主
请问各位大佬，mathematica在数值求解微分方程时需要根据求解的量，分成四个方程来求解：使用分段函数会报错：
在第一个参数 {<<3895220256 bytes>>} 中应该使用方程或者方程列表，而不是 \
Piecewise[<<3895218928 bytes>>
(*机架势能求解：每当换腿一次，机架的势能就要改一次*)
Subscript[m, w] = 0.730;
Subscript[F, T] = 0.8;
\[CapitalDelta]t = 0.001;
(*0-90°时前左腿*)
Subscript[y, G01] =
Subscript[r, 1]*Sin[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Sin[Subscript[\[Theta], ql2]] +
Subscript[r, 6]*Sin[Subscript[\[Theta], ql6]] +
Subscript[r, 9]*Sin[Subscript[\[Theta], ql9]];
(*机架势能*)
Subscript[EP,
01] = (Subscript[m, 0] + Subscript[m, w])*g*Subscript[y, G01];
(*用于求解各个杆件的势能，坐标系建在了曲柄转轴中心，距地面的势能为正的YG0-杆件质心高度位置*)
Subscript[EP,
I1] = -(2*Subscript[m, 1]*g*Subscript[y, G01] +
4*Subscript[m, 2]*g*Subscript[y, G01] +
4*Subscript[m, 3]*g*Subscript[y, G01] +
4*Subscript[m, 4]*g*Subscript[y, G01] +
4*Subscript[m, 5]*g*Subscript[y, G01]);
Subscript[x, G01] =
Subscript[r, 1]*Cos[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], ql2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], ql6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], ql9]];
Subscript[V, xG01] = D[Subscript[x, G01], t];
Subscript[V, yG01] = D[Subscript[y, G01], t];
(*机架的动能*)
Subscript[EK, 01] =
1/2*(Subscript[m, 0] + Subscript[m, w])*(Subscript[V, xG01]^2 +
Subscript[V, yG01]^2);
(*0-90总势能，动能*)
Subscript[EP, T1] =
Subscript[EP, I1] + Subscript[EP, 1] + Subscript[EP, 2] + Subscript[
EP, 3] + Subscript[EP, 4] - Subscript[EP, 01];
Subscript[EK, T1] =
Subscript[EK, 1] + Subscript[EK, 2] + Subscript[EK, 3] + Subscript[
EK, 4] + Subscript[EK, 01];
(*0-90度时曲柄等效力矩*)
Subscript[x, 1][t_] =
Subscript[r, 1]*Cos[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], ql2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], ql6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], ql9]];
\[CapitalDelta]xa[t_] :=
Subscript[x, 1][t + \[CapitalDelta]t] - Subscript[x, 1][t];
Subscript[\[Tau], 1] =
Subscript[F, T]*D[\[CapitalDelta]xa[t], Subscript[\[Theta], 1][t]];
a1 = D[Subscript[EP, T1], Subscript[\[Theta], 1][t]];
b1 = D[Subscript[EK, T1], Subscript[\[Theta], 1][t]];
c1 = D[D[Subscript[EK, T1], Subscript[\[Theta], 1]'[t]], t];
eq1 = c1 - b1 + a1 == Subscript[\[Tau], 1];
(*90-180度时*)
Subscript[y, G02] =
Subscript[r, 1]*Sin[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Sin[Subscript[\[Theta], qr2]] +
Subscript[r, 6]*Sin[Subscript[\[Theta], qr6]] +
Subscript[r, 9]*Sin[Subscript[\[Theta], qr9]];
Subscript[EP,
02] = (Subscript[m, 0] + Subscript[m, w])*g*Subscript[y, G02];
Subscript[EP,
I2] = -(2*Subscript[m, 1]*g*Subscript[y, G02] +
4*Subscript[m, 2]*g*Subscript[y, G02] +
4*Subscript[m, 3]*g*Subscript[y, G02] +
4*Subscript[m, 4]*g*Subscript[y, G02] +
4*Subscript[m, 5]*g*Subscript[y, G02]);
Subscript[x, G02] =
Subscript[r, 1]*Cos[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], qr2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], qr6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], qr9]];
Subscript[V, xG02] = D[Subscript[x, G02], t];
Subscript[V, yG02] = D[Subscript[y, G02], t];
Subscript[EK, 02] =
1/2*(Subscript[m, 0] + Subscript[m, w])*(Subscript[V, xG02]^2 +
Subscript[V, yG02]^2);
(*等效力矩*)
Subscript[x, 4][t_] =
Subscript[r, 1]*Cos[Subscript[\[Theta], 1][t]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], qr2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], qr6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], qr9]];
(*\[CapitalDelta]x2[t_]:=Subscript[x, \
4][t-\[CapitalDelta]t]-Subscript[x, 4][t];*)
\[CapitalDelta]xb[t_] :=
Subscript[x, 4][t + \[CapitalDelta]t] - Subscript[x, 4][t];
Subscript[\[Tau], 2] =
Subscript[F, T]*D[\[CapitalDelta]xb[t], Subscript[\[Theta], 1][t]];
(*90-180*)
Subscript[EP, T2] =
Subscript[EP, I2] + Subscript[EP, 1] + Subscript[EP, 2] + Subscript[
EP, 3] + Subscript[EP, 4] - Subscript[EP, 02];
Subscript[EK, T2] =
Subscript[EK, 1] + Subscript[EK, 2] + Subscript[EK, 3] + Subscript[
EK, 4] + Subscript[EK, 02];
a2 = D[Subscript[EP, T2], Subscript[\[Theta], 1][t]];
b2 = D[Subscript[EK, T2], Subscript[\[Theta], 1][t]];
c2 = D[D[Subscript[EK, T2], Subscript[\[Theta], 1]'[t]], t];
eq2 = c2 - b2 + a2 == Subscript[\[Tau], 2];
(*180-270du*)
Subscript[y, G031] =
Subscript[r, 1]*Sin[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Sin[Subscript[\[Theta], hl2]] +
Subscript[r, 6]*Sin[Subscript[\[Theta], hl6]] +
Subscript[r, 9]*Sin[Subscript[\[Theta], hl9]];
Subscript[EP,
03] = (Subscript[m, 0] + Subscript[m, w])*g*Subscript[y, G031];
Subscript[EP,
I3] = -(2*Subscript[m, 1]*g*Subscript[y, G031] +
4*Subscript[m, 2]*g*Subscript[y, G031] +
4*Subscript[m, 3]*g*Subscript[y, G031] +
4*Subscript[m, 4]*g*Subscript[y, G031] +
4*Subscript[m, 5]*g*Subscript[y, G031]);
Subscript[x, G031] =
Subscript[r, 1]*Cos[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], hl2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], hl6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], hl9]];
Subscript[V, xG03] = D[Subscript[x, G031], t];
Subscript[V, yG03] = D[Subscript[y, G031], t];
Subscript[EK, 03] =
1/2*(Subscript[m, 0] + Subscript[m, w])*(Subscript[V, xG03]^2 +
Subscript[V, yG03]^2);
Subscript[EP, T3] =
Subscript[EP, I3] + Subscript[EP, 1] + Subscript[EP, 2] + Subscript[
EP, 3] + Subscript[EP, 4] - Subscript[EP, 03];
Subscript[EK, T3] =
Subscript[EK, 1] + Subscript[EK, 2] + Subscript[EK, 3] + Subscript[
EK, 4] + Subscript[EK, 03];
Subscript[x, 2][t_] =
Subscript[r, 1]*Cos[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], hl2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], hl6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], hl9]];
\[CapitalDelta]xc[t_] :=
Subscript[x, 2][t + \[CapitalDelta]t] - Subscript[x, 2][t];
Subscript[\[Tau], 3] =
Subscript[F, T]*D[\[CapitalDelta]xc[t], Subscript[\[Theta], 1][t]];
a3 = D[Subscript[EP, T3], Subscript[\[Theta], 1][t]];
b3 = D[Subscript[EK, T3], Subscript[\[Theta], 1][t]];
c3 = D[D[Subscript[EK, T3], Subscript[\[Theta], 1]'[t]], t];
eq3 = c3 - b3 + a3 == Subscript[\[Tau], 3];
(*270-360度*)
Subscript[y, G04] =
Subscript[r, 1]*Sin[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Sin[Subscript[\[Theta], hr2]] +
Subscript[r, 6]*Sin[Subscript[\[Theta], hr6]] +
Subscript[r, 9]*Sin[Subscript[\[Theta], hr9]];
Subscript[EP,
04] = (Subscript[m, 0] + Subscript[m, w])*g*Subscript[y, G04];
Subscript[EP,
I4] = -(2*Subscript[m, 1]*g*Subscript[y, G04] +
4*Subscript[m, 2]*g*Subscript[y, G04] +
4*Subscript[m, 3]*g*Subscript[y, G04] +
4*Subscript[m, 4]*g*Subscript[y, G04] +
4*Subscript[m, 5]*g*Subscript[y, G04]);
Subscript[x, G04] =
Subscript[r, 1]*Cos[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], hr2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], hr6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], hr9]];
Subscript[V, xG04] = D[Subscript[x, G04], t];
Subscript[V, yG04] = D[Subscript[y, G04], t];
Subscript[EK, 04] =
1/2*(Subscript[m, 0] + Subscript[m, w])*(Subscript[V, xG04]^2 +
Subscript[V, yG04]^2);
Subscript[EP, T4] =
Subscript[EP, I4] + Subscript[EP, 1] + Subscript[EP, 2] + Subscript[
EP, 3] + Subscript[EP, 4] - Subscript[EP, 04];
Subscript[EK, T4] =
Subscript[EK, 1] + Subscript[EK, 2] + Subscript[EK, 3] + Subscript[
EK, 4] + Subscript[EK, 04];
Subscript[x, 3][t_] =
Subscript[r, 1]*Cos[Subscript[\[Theta], h1]] +
Subscript[r, 2]*Cos[Subscript[\[Theta], hr2]] +
Subscript[r, 6]*Cos[Subscript[\[Theta], hr6]] +
Subscript[r, 9]*Cos[Subscript[\[Theta], hr9]];
\[CapitalDelta]xd[t_] :=
Subscript[x, 3][t + \[CapitalDelta]t] - Subscript[x, 3][t];
Subscript[\[Tau], 4] =
Subscript[F, T]*D[\[CapitalDelta]xd[t], Subscript[\[Theta], 1][t]];
a4 = D[Subscript[EP, T4], Subscript[\[Theta], 1][t]];
b4 = D[Subscript[EK, T4], Subscript[\[Theta], 1][t]];
c4 = D[D[Subscript[EK, T4], Subscript[\[Theta], 1]'[t]], t];
eq4 = c4 - b4 + a4 == Subscript[\[Tau], 4];
(*动力学方程*)eqns =
Piecewise[{{eq1, -Pi/2 <= Subscript[\[Theta], 1][t] <= 0},
{eq2, -Pi <= Subscript[\[Theta], 1][t] < -Pi/2},
{eq3, -3*Pi/2 <= Subscript[\[Theta], 1][t] < -Pi},
{eq4, -2*Pi <= Subscript[\[Theta], 1][t] < -3*Pi/2}}];
(*eqns={Which[
-Pi/2<=Subscript[\[Theta], 1][t]<=0,eq1,
-Pi<=Subscript[\[Theta], 1][t]<-Pi/2,eq2,
-3*Pi/2<=Subscript[\[Theta], 1][t]<-Pi,eq3,
-2*Pi<=Subscript[\[Theta], 1][t]<-3*Pi/2,eq4
]};*)
solve1 =
NDSolve[{eqns,
WhenEvent[Abs[Subscript[\[Theta], 1][t]] >= 2*Pi - 1*Degree,
Subscript[\[Theta], 1][t] = 0], Subscript[\[Theta], 1][t] == 0,
Subscript[\[Theta], 1]'[0] == 0}, {Subscript[\[Theta], 1][t],
Subscript[\[Theta], 1]'[t], Subscript[\[Theta], 1]''[t]}, {t, 0,
3}, SolveDelayed -> True]
Plot[Evaluate[Subscript[\[Theta], 1][t] /. solve1], {t, 0, 3},
PlotRange -> All, AxesLabel -> {"t(s)", "曲柄角度(rad)"}]
Plot[Evaluate[Subscript[\[Theta], 1]'[t] /. solve1], {t, 0, 3},
PlotRange -> All, AxesLabel -> {"t(s)", "曲柄角速度(w)"}]
Plot[Evaluate[Subscript[\[Theta], 1]''[t] /. solve1], {t, 0, 3},
PlotRange -> All, AxesLabel -> {"t(s)", "曲柄角jia速度(w)"}]
data1 = Table[{t, Subscript[\[Theta], 1][t] /. solve1,
Subscript[\[Theta], 1]'[t] /. solve1,
Subscript[\[Theta], 1]''[t] /. solve1}, {t, 0, 3, 0.001}]
Export["output.txt", data1, "Table"]
结果总是报错
2025年01月02日 14点01分 1
吧务
level 15
xzcyr
1. NDSolve至少现在还不支持你想象中的那个语法（一个Piecewise里塞几个等式这种），请把表达式塞到Piecewise里再构造方程。这部分最简单的改法大概是：
eqns = Piecewise[{{Subtract @@ eq1, -(π/2) <=
Subscript[θ, 1][t] <= 0}, {Subtract @@ eq2, -π <=
Subscript[θ, 1][t] < -(π/2)}, {Subtract @@
eq3, -((3 π)/2) <=
Subscript[θ, 1][t] < -π}, {Subtract @@
eq4, -2 π <=
Subscript[θ, 1][t] < -((3 π)/2)}}] == 0;
2. 初始条件有低级错误。
3. 若干个 r 和 m 的下标变量定义没给。
全改过来就能算了。
2025年01月04日 12点01分 2
1