求教ode45求解lorenz方程

jlliang

用ode45求解lorenz方程，只是将时间设为Tspan=[0  5],算了半天都没算出来，是什么原因？
记得以前算还挺快的，现在为何变慢?我的代码如下：

[t,y]=ode45(@lorenz,[0 5],[15.34;13.68;37.91]);

function dy=lorenz(t,y0)
%a=10,b=8/3;c=28;
%dx/dt=a(y-z),
%dy/dt=cz-xz-y,
%dz/dt=xy-bz
a=10;b=8/3;c=28;
dy=zeros(3,1);
dy(1,1)=a*(y0(2)-y0(3));
dy(2,1)=c*y0(3)-y0(1)*y0(3)-y0(2);
dy(3,1)=y0(1)*y0(2)-b*y0(3);

[ 本帖最后由 jlliang 于 2009-10-17 13:52 编辑 ]

无水1324

测试了一下，确实比较奇怪！
怎么会出现这种问题呢

jlliang

方程找到了，不好意思
方程应该是
%a=10,b=8/3;c=28;
%dx/dt=a(y-x),
%dy/dt=cx-xz-y,
%dz/dt=xy-bz

octopussheng

楼主可以试试分析分析，为什么方程仅仅变了一个量，计算速度就差这么多，与大家共享哈。

无水1324

确实，可以从这个角度考虑一下。

jlliang

我觉得与方程的特性有关，ode45不适合刚性问题，且ode45是用变步长的的-4阶R-K法，绝对误差为e-3,相对误差为e-6（默认值），数值解无法满足精度要求，便步长消耗了时间。这是个人想法，需要验证。

octopussheng

其实ode45也不一定就是变步长，它的计算步长也可以设定为定步长计算的。

jlliang

ode45的定步长应该是这样设置的吧[0:0.01:100],  步长是t＝0.01吧，如果是这样的，我有点疑问，计算出的是该点处的数值解，但这样的步长会满足它的默认精度吗，如果不满足，该点处的解又是怎么得到的呢？是减少步长使其满足精度，并能求出t＝0.01处的数值解，还是..........?
实际中，在我们设定精度条件下，仍能求出指定时刻的数值解，这样的话，定步长t＝？

octopussheng

在ode45计算的时候，精度仍是可以设置的，如果不满足，会有警告信息。
如果不满足的话，只有减少时间步长了。

你最后一句话说的意思我不大明白，请解释一下。

jlliang

比如设定精度options＝odeset（'Reltol',1e-10,'Abstol',1e-10），Tspan=[0:0.1:100]
我要的是指定时刻 t=n*0.1时刻的数值解，那么定步长是h=0.1，还是满足精度下的步长（有两种情况h＝0.1，一种比0.1小），如果是比0.1小，那么我们是不能设定定步长计算的
这样就不是真正的定步长计算了

free518

没计算过;学习了!!!!!!!!!!!!!!!!

octopussheng

如果真实步长比0.1小，那么计算就不能满足其仿真误差要求，那么只有降低步长，才能满足计算要求了。