2

質量ばねと二重振り子を備えた (やや奇妙な) システムの運動方程式をシミュレートしています。このシステムには、質量行列と関数 f(x) があり、ode45 を呼び出して解決します。

M*x' = f(x,t);

5 つの状態変数 q= [ QDot, phi, phiDot, r, rDot]'; があります。(QDot は現在の QDot に依存するものがないため削除されました。) ここで、いくつかの力を計算するために、ode45 が各積分ステップで計算する rDotDot の計算値も保存したいと思いますが、ode45 はこれを返しません。私は少し調べましたが、私が見つけた唯一の2つの解決策は、a）これを3次問題に変え、状態ベクトルにphiDotDotとrDotDotを追加することです。これは可能な限り回避したいと思います。これはすでに非線形であり、これにより事態がさら​​に悪化し、計算時間が大幅に増加するからです。

b)ここで説明されているように、状態を拡張して関数を直接計算します。ただし、例では、質量行列にゼロの行を追加すると彼は言います。そうしないと導関数を積分し、ある点で評価するだけでなく、質量行列を特異にするので、これは理にかなっています。私にはうまくいかないようです...

これは基本的なことのように思えます (状態ベクトルの微分値を取得するため)。私が考えていない明らかなことはありますか? （または、それほど明白ではないものでも構いません....）

ああ、そしてグローバル変数はそれほど素晴らしいものではありません。なぜなら、ode45 は f() 関数をステップの調整中に何度も呼び出すため、グローバル変数のサイズと返される状態ベクトル q がまったく一致しないからです。

誰かがそれを必要とする場合に備えて、コードは次のとおりです。

%(Initialization of parameters are above this line)
   options = odeset('Mass',@massMatrix);
   [T,q] = ode45(@f, tspan,q0,options);

function dqdt = f(t,q,p)
    % q = [qDot phi phiDot r rDot]';

    dqdt = zeros(size(q));

    dqdt(1) = -R/L*q(1) - kb/L*q(3) +vs/L;
    dqdt(2) = q(3);
    dqdt(3) = kt*q(1) + mp*sin(q(2))*lp*g;
    dqdt(4) = q(5);
    dqdt(5) = mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g;
end

function M = massMatrix(~,q)
    M = [
        1 0 0 0 0;
        0 1 0 0 0;
        0 0 mp*lp^2 0 -mp*lp*sin(q(2));
        0 0 0 1 0;
        0 0 mp*lp*sin(q(2)) 0 (mb+mp)
        ];
end
4

2 に答える 2

2

ネストされた関数を使用しているため、スコープ規則を使用してグローバル変数の動作を模倣できます。

この目的のために出力関数を (ab) 使用するのが最も簡単です。

function solveODE

    % ....        
    %(Initialization of parameters are above this line)

    % initialize "global" variable
    rDotDot = [];

    % Specify output function 
    options = odeset(...
        'Mass', @massMatrix,...
        'OutputFcn', @outputFcn);

    % solve ODE
    [T,q] = ode45(@f, tspan,q0,options);

    % show the rDotDots    
    rDotDot



    % derivative 
    function dqdt = f(~,q)

        % q = [qDot phi phiDot r rDot]';

        dqdt = [...
            -R/L*q(1) - kb/L*q(3) + vs/L
            q(3)
            kt*q(1) + mp*sin(q(2))*lp*g
            q(5)
            mp*lp*cos(q(2))*q(3)^2 - ks*q(4) - (mb+mp)*g
        ];

    end % q-dot function 

    % mass matrix
    function M = massMatrix(~,q)
        M = [
            1 0 0 0 0;
            0 1 0 0 0;
            0 0 mp*lp^2 0 -mp*lp*sin(q(2));
            0 0 0 1 0;
            0 0 mp*lp*sin(q(2)) 0 (mb+mp)
         ];
    end % mass matrix function


    % the output function collects values for rDotDot at the initial step 
    % and each sucessful step
    function status = outputFcn(t,q,flag)

        status = 0;

        % at initialization, and after each succesful step
        if isempty(flag) || strcmp(flag, 'init')
            deriv = f(t,q);
            rDotDot(end+1) = deriv(end);
        end

    end % output function 

end

出力関数は、最初のステップと成功したすべてのステップで導関数のみを計算するため、基本的に Adrian Ratnapala が提案したのと同じことを行っています。の各出力で導関数を再実行しますode45。それはさらにエレガントだと思います（エイドリアンの場合は+1）。

出力関数アプローチには、統合の実行中に値をプロットできるという利点がありますrDotDot(! を忘れないでくださいdrawnow)。これは、長時間実行される統合に非常に役立ちます。

于 2013-05-14T09:34:20.560 に答える
2

最も簡単な解決策は、 によって返された各値に対して関数を再実行することode45です。

難しい解決策は、DotDots を別の場所 (事前に割り当てられたマトリックスまたは外部ファイル) に記録することです。問題は、ode45 が奇妙な場所で密かに評価を行うと、不要なデータ ポイントが作成される可能性があることです。

于 2013-05-14T08:31:11.380 に答える