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'))

必要に応じて、2 つの方法のいずれかで警告メッセージを回避することができます。

1.警告をオフにします (後で再度有効にしてください)。次のコードでそれを行うことができます。

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

関数を使用して、waring の ID を取得できlastwarnます。

2.もう 1 つのオプションは、統合の境界を変更し、警告を完全に回避することです。たとえば、次のようになります。

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

これは許容できる場合と許容できない場合がありますがintegral、結果が収束しないため、正しい解が返されません。

于 2014-02-01T17:16:31.057 に答える