3

MATLAB を使用して、次のような方程式を解こうとしています。

B = alpha*Y0*sqrt(epsilon)/(pi*ln(b/a)*sqrt(epsilon_t))*0 から pi までの (2*sinint(k0*sqrt(epsilon*(a^2+b) ^2-2abcos(シータ))-sinint(2*k0*sqrt(イプシロン)*a*sin(シータ/2))-sinint(2*k0*sqrt(イプシロン)*b*sin(シータ/2)) ) シータに関して

イプシロンは未知数です。

と を使用して積分に未知数が埋め込まれた方程式をシンボリックに解く方法を知っていますが、シンボリック積分器を使用するint()とsolve()、int()この複雑な方程式には時間がかかりすぎます。quad()、quadl()およびを使用しようとするとquadgk()、未知数が積分にどのように埋め込まれているかを処理するのに苦労します。

4

1 に答える 1

2

この種のことは、本当に急速に複雑になります。単一のインライン式ですべてを実行することは可能ですが、読みやすくするためだけであれば、複数のネストされた関数に分割することをお勧めします。

読みやすさが重要である理由の最も良い例: 投稿した数式にブラケットの問題があります。閉じ括弧が足りないので、方程式が数学表記でどのように見えるか完全にはわかりません:)

とにかく、これがあなたが意図したバージョンでそれを行う1つの方法です。

function test

    % some random values for testing    

    Y0 = rand;
    b  = rand;
    a  = rand;
    k0 = rand;
    alpha     = rand;
    epsilon_t = rand;

    % D is your B
    D = -0.015;

    % define SIMPLE anonymous function 
    Bb = @(ep) F(ep).*main_integral(ep) - D;

    % aaaand...solve it!
    sol = fsolve(Bb, 1)

    % The anonymous function above is only simple, because of these: 

    % the main integral    
    function val = main_integral(epsilon)

        % we need for loop through epsilon, due to how quad(gk) solves things
        val = zeros(size(epsilon));
        for ii = 1:numel(epsilon)

            ep = epsilon(ii);

            % NOTE how the sinint's all have a different function as argument:
            val(ii) = quadgk(@(th)...
                2*sinint(A(ep,th)) - sinint(B(ep,th)) - sinint(C(ep,th)), ...
                0, pi);            
        end

    end

    % factor in front of integral
    function f = F(epsilon)        
        f = alpha*Y0*sqrt(epsilon)./(pi*log(b/a)*sqrt(epsilon_t)); end

    % first sinint argument
    function val = A(epsilon, theta)
        val = k0*sqrt(epsilon*(a^2+b^2-2*a*b*cos(theta))); end

    % second sinint argument
    function val = B(epsilon, theta)
        val = 2*k0*sqrt(epsilon)*a*sin(theta/2); end

    % third sinint argument
    function val = C(epsilon, theta)
        val = 2*k0*sqrt(epsilon)*b*sin(theta/2); end    

end

上記の解決策は依然として非常に遅くなりますが、これほど複雑な積分ではかなり普通のことだと思います。

sinint速度の低下のほとんどは非組み込み関数を使用した for ループによるものであるため、独自の実装はあまり役に立たないと思います...速度が必要な場合は、独自の Gauss を使用した MEX 実装を使用します-クロンロッド適応直交ルーチン。

于 2012-12-06T12:54:01.973 に答える