Как использовать MATLAB для численного решения уравнения с неизвестным, встроенным в интеграл?

Я пытался использовать 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(), у меня возникают проблемы с пониманием того, как неизвестное встраивается в интеграл.


matlab equation integral
person user1880273    schedule 05.12.2012    source источник
comment
mathworks.com/matlabcentral/fileexchange/   -  person 0x90    schedule 06.12.2012
comment
Все ли остальные константы определены выше? (например, альфа, Y0 и т. д.)   -  person nicktruesdale    schedule 06.12.2012
comment
Эти константы являются просто произвольными константами для целей расчета, только эпсилон является неизвестным. тета - это просто переменная интегрирования (например, dx, dy, dz и т. д.)   -  person user1880273    schedule 06.12.2012

Ответы (1)


arrow_upward
2
arrow_downward

Такие вещи очень быстро усложняются. Хотя все это можно сделать в одном встроенном уравнении, я бы посоветовал вам разбить его на несколько вложенных функций, хотя бы для удобства чтения.

Лучший пример того, почему удобочитаемость важна: у вас есть проблема с скобками в опубликованном вами уравнении; закрывающих скобок недостаточно, поэтому я не могу быть полностью уверен, как выглядит уравнение в математической записи :)

Во всяком случае, вот один из способов сделать это с версией, которую я -- думаю -- вы имели в виду:

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 с не встроенными функциями... Если вам нужна скорость, я бы пошел на реализацию MEX с вашим собственным Адаптивная квадратурная процедура Гаусса-Кронрода.

person Rody Oldenhuis    schedule 06.12.2012