for循环可能出了问题
mathematica吧
全部回复
仅看楼主
level 2
- 楼主
大佬们,想用mathematica解一个二阶微分方程(微振动),画出x(t)和x'(t)的图像;但是t如果太大,图像会“飞”出去,所以想用for循环,利用能量守恒为判据,来确定最大的t,但不知道为什么会报错...第一次用mathematica...
x[0] = 0
x'[0] = 0.4
C = 0.5 x'[0]^2 + 1/3 x[0]^3 + 0.5 x[0]^2
For[i = 0, 0.5 x'[i]^2 + 1/3 x[i]^3 + 0.5 x[i]^2 <= C, i + 0.01,
NDSolve[{x''[t] + x[t] + x[t]^2 == 0, x[0] == 0, x'[0] == 0.6},
x, {t, 0, i}]]
ParametricPlot[Evaluate[{x[t], x'[t]} /. %], {t, 0, i}]
2022年03月31日 04点03分 1
吧务
level 10
一般WhenEvent StopIntegration来终止计算
2022年03月31日 08点03分 2
此外,你的一些赋值行为会带来不良后果。C的赋值是不可行的因为它是内部符号,对导数的赋值不易清除。
2022年03月31日 08点03分
-
非常感谢,我再好好研究研究mathematica
2022年03月31日 16点03分
LZ这个根本轮不到WhenEvent,t=9.86多的时候就撞上奇点自己停了。
2022年04月02日 02点04分
吧务
level 15
你这思路就离谱。全算完再找符合条件的t不就完了?
xsol = NDSolveValue[{x''[t] + x[t] + x[t]^2 == 0, x[0] == 0, x'[0] == 6/10},
x, {t, 0, 10}]
1/2 x'[t]^2 + 1/3 x[t]^3 + 1/2 x[t]^2 /. x -> xsol /. t -> 0
Plot[1/2 x'[t]^2 + 1/3 x[t]^3 + 1/2 x[t]^2 /. x -> xsol, {t, 0, 9.3}, PlotRange -> All]
(* 尽管图像是这个样子,但我并不觉得这是因为能量不守恒了,因为奇点附近的计算精度肯定会受影响。姑且试着算下能量超出1%的位置: *)
FindRoot[1.01 0.18 == 1/2 x'[t]^2 + 1/3 x[t]^3 + 1/2 x[t]^2 /. x -> xsol, {t, 8}]
(* {t -> 9.4571} *)
2022年04月01日 16点04分 3
-
嗯...我们老师也说是因为计算精度的问题,会导致图像“飞”出去。我再想想,谢谢大佬指点
2022年04月02日 05点04分
1