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(theta))-sinint(2*k0*sqrt(epsilon)*a*sin(theta/2))-sinint(2*k0*sqrt(epsilon)*b*sin(theta/2)) ) 关于 theta

epsilon是未知数。

我知道如何使用int()and符号地求解嵌入在积分中的未知方程solve()，但是对于如此复杂的方程，使用符号积分器int()需要太长时间。当我尝试使用quad(), quadl()andquadgk()时，我无法处理未知数如何嵌入积分中。

4

1 回答 1

2

这种事情变得非常复杂。尽管可以在一个内联方程中完成所有操作，但我建议您将其拆分为多个嵌套函数，如果只是为了便于阅读。

为什么可读性很重要的最好例子：你在你发布的方程式中有一个括号问题；没有足够的右括号，所以我不能完全确定等式在数学符号中的样子:)

无论如何，这是使用我认为的版本的一种方法 - 你的意思是：

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 -Kronrod 自适应正交例程。

于 2012-12-06T12:54:01.973 回答