在看MATLAB中ode45的应用

雨人

各位高手：在应用matlab中ODE45求解二阶微分方程时，如何输出状态方程本身随时间的时间历程。比如方程是二阶常微分方程，应用ODE45可以求出X和X的一次导数,如何获得X的二次导数呢？请高手们指导？

雨人

解析解：
time=linspace(0,4,2000);
x=dsolve('5*D2x+2000*x+0.5*5*9.81*Dx=100*sin(30*t)','Dx(0)=0.1','x(0)=0.1');
xdot=diff(x);
x2dot=diff(xdot);

数值方法：
function solve_eq_examp
global k acc;
k=0;    % 该计数器输出状态方程被调用的次数
acc=[]; %输出状态方程的时间历程
[t x]=ode45(@Slide_equ,[0:0.01:4],[0.1 0.1]);
子函数
function xdot=Slide_equ(t,x)
global k acc;
k=k+1
xdot=[x(2)
      -2000/5*x(1)-0.5*9.8*sign(x(2))+100/5*sin(30*t)];
  acc=[acc;[t xdot']]; %这里输出状态方程的时间历程
  

雨人

数值方法中该计数器最后为3385，表明状态方程被调用了3385次。在应用全局变量返回状态方程每次被调用的是的结果后，得到的x二阶导数时间历程和解析解差距较大？
全局变量acc返回的矩阵大小为3385向，要不ODE放回的[t x] 大很多？？

雨人

请问高手，有什么好的解决方法？不胜感激？

雨人

怎么一直没有人关注啊！:@L

gghhjj

雨人 发表于 2011-3-26 12:33

                               
登录/注册后可看大图

数值方法中该计数器最后为3385，表明状态方程被调用了3385次。在应用全局变量返回状态方程每次被调用的是的 ...

说明一下你是如何求二阶导数的？怀疑可能是数值求导方面的问题

最好把代码给全，并附上结果图

tudexin

!

meiyongyuandeze

请问X何X一阶导数和解析值相差大吗？ODE45也不是所有的微分方程都是用！

zhao522

回复 8 # meiyongyuandeze 的帖子

ode45是用来求解二阶微分方程，具体的实现是通过变量代换，将二价的转换为一阶的进行求解的，方程的输出只有X和X的倒数，你如果想得到X的二阶倒数，你可以根据方程本身，将X的二阶倒放到一边，其他的放到另一边，直接代入数值就可以得到X的二阶倒，也不知道这样可不可行。

雨人

E:\matlab\Steerability\单自由度方程的不同解法、解析数值解对比.bmp

感谢大家的回复，以上代码是完整的，可以使用，这里我把结果贴上，以供大家讨论。

雨人

E:\matlab\Steerability\单自由度方程的不同解法\解析数值解对比.bmp

我上传不了图片！:@L

        对于之前发的帖子，现在有点结果了。做这个简单的方程是为了验证一下方法的可行性。由于ODE45在积分时，多次调用了被积函数，直接使用全局变量返回被积函数的每次迭代值是不准确的。图片上的加速度项可以看出误差的波动项。
        之前之所以想使用二阶微分项的原因是，由于动力系统各个自由度之间存在着惯性耦合，(同一个微分方程中包含着其它自由度上微分方程变量的二阶微分项)，所以进行数值求解时遇到了困难，如果使用SIMULINK去左，出现代数环，同样很难解决。现在的决定从方程本身着手，手动解除微分方程组的二阶导数之间的耦合项，一遍求解方便。
再次谢谢大家的支持

雨人

回复 6 # gghhjj 的帖子

图片上传出错了，程序是完整的，感谢参与讨论

meiyongyuandeze

关注，希望能分享的进展！

雨人

回复 13 # meiyongyuandeze 的帖子

个人感觉使用ODE系列函数进行数值求解，也就是对状态方程的积分，获得各个状态变量的时间序列。由于ODE系列方程只能返回各个状态变量积分后的时间序列，我的问题有什么方法可以把各个时间点对应的状态也输出出来。这个问题我一直没有解决。考虑过全局变量，但是结果不理想。也使用过用积分后的状态变量重新计算带入状态方程，计算各个时间点的状态。感觉有点麻烦。请问各位有什么好的方法。

雨人

                               
登录/注册后可看大图