CODEWITHSUNDEEP
matlabequationintegral
user1880273
user1880273

Reputation: 31

How to use MATLAB to numerically solve equation with unknown embedded in integral?

I've been trying to use MATLAB to solve equations like this:

B = alpha*Y0*sqrt(epsilon)/(pi*ln(b/a)*sqrt(epsilon_t))*integral from 0 to pi of (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))) with regard to theta

Where epsilon is the unknown.

I know how to symbolically solve equations with unknown embedded in an integral by using int() and solve(), but using the symbolic integrator int() takes too long for equations this complicated. When I try to use quad(), quadl() and quadgk(), I have trouble dealing with how the unknown is embedded in the integral.

Upvotes: 3

Views: 4178

Answers (1)

Rody Oldenhuis
Rody Oldenhuis

Reputation: 38032

This sort of thing gets complicated real fast. Although it is possible to do it all in a single inline equation, I would advise you to split it up into multiple nested functions, if only for readability.

The best example of why readability is important: you have a bracketing problem in the eqution you posted; there's not enough closing brackets, so I can't be entirely sure what the equation looks like in mathematical notation :)

Anyway, here's one way to do it with the version I --think-- you meant: 

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

The solution above will still be quite slow, but I think that's pretty normal for integrals this complicated.

I don't think implementing your own sinint will help much, as most of the speed loss is due to the for loops with non-builtin functions...If it's speed you want, I'd go for a MEX implementation with your own Gauss-Kronrod adaptive quadrature routine.

Upvotes: 2

Related Questions