maple或maplesim似乎不能求解带有事件的常微分方程？一道小题

TBE_Legend

maple或maplesim似乎不能求解带有事件的常微分方程？

比如： mathematica自带的这个例子， 从静摩擦到动摩擦的模拟：

我个人能力有限，似乎maple或maplesim解不了这个方程。

feiyuzhen

前三个模型的求解代码如下。

1. 
   
2. restart;
   
3. m, l := 1., 1.;
   
4. beltv := proc (t) options operator, arrow; .1 end proc;
   
5. spring := proc (x) options operator, arrow; 1000.*(l-x) end proc;
   
6. sys1 := [m*(diff(x(t), t, t)) = spring(x(t))+friction1(diff(x(t), t)), x(0) = 1, (D(x))(0) = 0];
   
7. viscous := proc (v) options operator, arrow; (-1)*30.*(v-beltv(t)) end proc;
   
8. friction1 := proc (v) options operator, arrow; viscous(v) end proc;
   
9. friction1(v);
   
10. sys1;
    
11. sol1 := dsolve(sys1, numeric, output = listprocedure, range = 0 .. 2);
    
12. plots[odeplot](sol1, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2, refine = 2);
    
13. plots[odeplot](sol1, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, refine = 2);
    
14. sys2 := [m*(diff(x(t), t, t)) = spring(x(t))+friction2(diff(x(t), t)), x(0) = 1, (D(x))(0) = 0];
    
15. coulomb := proc (v) options operator, arrow; (-1)*25.*signum(v-beltv(t)) end proc;
    
16. friction2 := proc (v) options operator, arrow; viscous(v)+coulomb(v) end proc;
    
17. friction2(v);
    
18. sys2;
    
19. sol2 := dsolve(sys2, numeric, method = classical[abmoulton], stepsize = 0.4500000e-4, output = listprocedure);
    
20. plots[odeplot](sol2, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2);
    
21. plots[odeplot](sol2, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, numpoints = 500);
    
22. sys3 := [m*(diff(x(t), t, t)) = spring(x(t))+friction3(diff(x(t), t)), x(0) = 1, (D(x))(0) = 0];
    
23. stribeck := proc (v) options operator, arrow; (-1)*.3*signum(v)*exp(-2*abs(v)) end proc;
    
24. friction3 := proc (v) options operator, arrow; viscous(v)+coulomb(v)+stribeck(v) end proc;
    
25. friction3(v);
    
26. sys3;
    
27. sol3 := dsolve(sys3, numeric, method = classical[abmoulton], stepsize = 0.4500000e-4, output = listprocedure);
    
28. plots[odeplot](sol3, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2);
    
29. plots[odeplot](sol3, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, numpoints = 500);
    
30. plots[display]({plots[odeplot](sol1, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2, refine = 2), plots[odeplot](sol2, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2, color = green), plots[odeplot](sol3, [t, x(t)], 0 .. 1, gridlines = true, thickness = 2, color = blue)}, view = [0 .. 1, 1 .. 1.03]);
    
31. plots[display]({plots[odeplot](sol1, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, refine = 2), plots[odeplot](sol2, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, numpoints = 500, color = green), plots[odeplot](sol3, [[t, diff(x(t), t)], [t, beltv(t)]], 0 .. 1, gridlines = true, thickness = 2, numpoints = 500, color = blue)}, view = [0 .. 1, -0.2e-1 .. .13]);

复制代码

feiyuzhen

个人感觉能求，不过对带事件的常微分方程不熟，
第一个方程好算，中间两个也能算，不过初始时间段似乎有点问题。

TBE_Legend

feiyuzhen 发表于 2013-3-29 08:58
个人感觉能求，不过对带事件的常微分方程不熟，
第一个方程好算，中间两个也能算，不过初始时间段似乎有点 ...

我试了下 comsol 和 simulationx 都可以解。

simx 用了 typedesigner。

TBE_Legend

mathematica 的确很好用， 但有bug， 你用用就知道了， 不信你把 belt 的速度改为与时间相关的。

comsol 也很好，以后可以和pde一起做。

TBE_Legend

以上的结果都是 对应于 mathematica 自带例子中 第一个 friction model，就是最简单的那个。

不过希望可以看看maplesim怎么做，正在下载 maplesim v6 中， 有空在研究吧

TBE_Legend

真不错，多谢老兄啊，看来maole是可以的，我回去仔细研究下先。