求教怎么在数值求解微分方程组的过程中添加导数的初始条件?
mathematica吧
全部回复
仅看楼主
level 2
tim1104a 楼主
我在求解微分方程组的过程中试图给函数的导数设置初始条件,使用了如y'[0]==0、D[y[t], t] /. t -> 0 == 0这样的方法,但总是报错:应该使用方程或者方程列表,而不是 y'[True],请问这是怎么回事呢?
下面是代码:
Clear[x1, x2, y1, y2, z1];
deqns = {
m1*x1''[t] ==
m2*x2''[t]*(
l - Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])/(
x1[t] - x2[t])*((\[Delta] Cos[\[Omega] t] - x1[t])/
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2] + (
x2[t] - x1[t])/(
l - Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])),
m1*y1''[t] ==
m2*x2''[t]*(
l - Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])/(
x1[t] - x2[t])*((\[Delta] Sin[\[Omega] t] - y1[t])/
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2] + (
y2[t] - y1[t])/(
l - Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])),
m1*z1''[t] ==
m2*x2''[t]*((-(l -
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])*z1[t])/(
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2]*(x1[t] -
x2[t])) +
Sqrt[(l -
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])^2 - (x2[
t] - x1[t])^2 - (y2[t] - y1[t])^2]/(x1[t] - x2[t])) + m1*g,
m2*y2''[t] == m2*x2''[t]*(y1[t] - y2[t])/(x1[t] - x2[t]),
m2*Dt[Sqrt[(l -
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])^2 - (x2[
t] - x1[t])^2 - (y2[t] - y1[t])^2] + z1[t], {t, 2}] ==
m2*x2''[t]*-Sqrt[(l -
Sqrt[(\[Delta] Sin[\[Omega] t] - y1[t])^2 +
z1[t]^2 + (\[Delta] Cos[\[Omega] t] - x1[t])^2])^2 - (x2[
t] - x1[t])^2 - (y2[t] - y1[t])^2]/(x1[t] - x2[t]) + m2*g
};
ics = {
x1[0] == 1/4 l + \[Delta], y1[0] == 0, z1[0] == Sqrt[3]/4 l,
D[x1[t], t] /. t -> 0 == 0, D[y1[t], t] /. t -> 0 == v0,
D[z1[t], t] /. t -> 0 == 0,
x2[0] == (1 - Sqrt[3])/4 l + \[Delta], y2[0] == 0(*,z2[0]\[Equal](
Sqrt[3]+1)/4l
x2'[0]\[Equal]0,y2'[0]\[Equal]0,z2'[0]\[Equal]0*)
};
params = {m1 -> 0.023, m2 -> 0.023, l -> 0.6,
g -> 9.81, \[Omega] -> 20, \[Delta] -> 0.05, v0 -> 1};
res = NDSolve[{deqns~Join~ics} /. params, {x1, y1, x2, y2, z1}, {t, 0,
15}, Method -> {"EquationSimplification" -> "Residual"}]
2020年02月06日 15点02分 1
level 2
tim1104a 楼主
已解决,可能是我之前写过y'[0]=0这样的语句,重启一下就没问题了
2020年02月07日 02点02分 2
y'[0]的FullForm是Derivative[1][y][0],所以只需Clear[Derivative]即可清除其定义
2020年02月07日 04点02分
1