0

当我尝试使用这个 Matlab 代码时，它进入了一个无限循环。我正在尝试在内部执行集成ode45：

clear
clc
options = odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);

plot(T,Y(:,1),'+',T,Y(:,2),'*',T,Y(:,3),'.')


function dy = rigid(t,y)
dy = zeros(3,1);    % a column vector

dy(1) = y(2) ;

dy(2) = -y(1) * y(3);
fun = @(t) exp(-t.^2).*log(t).^2+y(1);
q = integral(fun,0,Inf);
dy(3) = y(2) * y(3) + q;
4

1 回答 1

0

没有“无限循环”。您的功能只需要很长时间才能集成。尝试设置tspan为[0 1e-7]. 它似乎是一个高频振荡，但我不知道你的方程是否正确（这是一个数学问题而不是编程问题）。这样的系统很难准确集成（ode15可能是更好的选择），更不用说快速集成了。

您也没有费心提及调用integral生成警告消息的重要事实：

Warning: Minimum step size reached near x = 1.75484e+22. There may be a
singularity, or the tolerances may be too tight for this problem. 
> In funfun/private/integralCalc>checkSpacing at 457
  In funfun/private/integralCalc>iterateScalarValued at 320
  In funfun/private/integralCalc>vadapt at 133
  In funfun/private/integralCalc at 84
  In integral at 88
  In rtest1>rigid at 17
  In ode15s at 580
  In rtest1 at 5

在每次迭代中打印出警告消息会大大减慢集成速度。这个警告有充分的理由。您确实意识到您使用integralfrom 0to评估的功能Inf等同于以下内容，对吗？

sqrt(pi)*((eulergamma + log(4))^2/8 + pi^2/16) + Inf*y(1)

或eulergamma在哪里。你的被积函数不收敛。psi(1)double(sym('eulergamma'))

如果您愿意，可以尝试通过以下两种方式之一来避免警告消息。

1.关闭警告（确保之后重新启用它）。您可以使用以下代码执行此操作：

...
warning('OFF','MATLAB:integral:MinStepSize');
[T,Y] = ode45(@rigid,[0 12],[0 1 1],options);
warning('ON','MATLAB:integral:MinStepSize');
...

lastwarn您可以通过该函数获取waring的ID 。

2.另一种选择可能是更改您的集成范围并完全避免警告，例如：

...
q = integral(fun,0,1e20);
...

这可能是可接受的，也可能是不可接受的，但integral也不会返回正确的解决方案，因为结果不会收敛。

于 2014-02-01T17:16:31.057 回答