MATLAB: How can i calculate double integral

integralsyms

Good afternon! I have the next problem: I try to calculate double integral, but I receive next error message: 
Error using integral2Calc>integral2t/tensor (line 231)
Input function must return 'double' or 'single' values. Found 'symfun'.
Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);
Error in integral2Calc (line 9)
    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in Integra (line 21)
    q = integral2(innerI, 0, 2*pi, lim2, pi/2);
This is full code. Can some one help me?
function innerI = Integra(f, num)
    syms theta phi
    global theta phi w rho bound kb kOm epsROm a dd Rx THx FIx;
    rho = 0.04; %радиус сферы
    bound = 0.045; %радиус границы рассеивателя
    dd = 10; %расстояние от границы до возмущения
    a = 0.1; %радиус обзора возмущения
    Rx = 5; THx = pi/4; FIx = 3*pi/4;
    epsROm = 4; muROm = 1; sigmaOm = 0.01; %Параметры среды
    epsRB = 6;  muRB = 1;  sigmaRB = 1; %Параметры возмущения
    w = 2*pi*f; %Круговая частота
    kb = wave_number(muRB, epsRB, sigmaRB);
    kOm = wave_number(muROm, epsROm, sigmaOm);
    nmax = num;
    lim2 = pi/2 - atan(a/dd);
    innerI = @(theta, phi) abs(Escr(nmax) * gradFk()) - abs(Esct(nmax) * gradFk()) * ...
        (-dd*cos(theta)/(sin(theta) * sin(theta)));
    q = integral2(innerI, 0, 2*pi, lim2, pi/2);
    % res = int(innerI, phi, [0, 2*pi]);
%     re_2 = int(res, theta, [lim2, pi/2]);
%     disp(re_2);
end
function k = wave_number(mu, eps, sigma)
    global w;
    mu0 = 4*pi*10^(-7);
    eps0 = 8.8541878128*10^(-12);
    k = w * sqrt(mu * mu0 * eps0 * (eps + 1i *(sigma /(w * eps0))));
end
function Escrr = Escr(nmax)
    global phi bound kOm;
    k_rho = kOm*bound;
    sum = 0;
    for n = 1:nmax
        [Ben, ~] = get_coeffs(n);
        sum = sum + n*(n+1) * Ben * Dzeta(n, k_rho) * lejandr(n);
    end
    Escrr = sum * cos(phi)/(k_rho)^2;
end
function Esctt = Esct(nmax)
    global theta phi bound kOm;
    k_rho = kOm*bound;
    sum = 0;
    for n = 1:nmax
        [Ben, Bmn] = get_coeffs(n);
        sum = sum + Ben * Dzeta_1(n, k_rho)*lejandr_1(n)*sin(theta) - ...
            1i * Bmn * Dzeta(n,k_rho) * lejandr(n)/sin(theta);
    end
    Esctt = (-cos(phi)/(k_rho))*sum;
end
function Fik = gradFk()
    global theta phi rho kOm Rx THx FIx;
    syms d Fik
    d = dist(THx, FIx, rho, Rx);
    denominator = 8 * d(theta, phi);
    fact1 = 1i * kOm;
    fact2 = -2 * besselh(1, kOm*d(theta, phi));
    fact3 = rho-Rx*(cos(theta-THx)*sin(FIx)*sin(phi)+cos(FIx)*cos(phi));
    Fik = fact1 * fact2 * fact3 / denominator;
end
function [Ben, Bmn] = get_coeffs(n)
    global w rho kb kOm epsROm;
    c = 3*10^8;
    q = w*rho*sqrt(epsROm)/c;
    N = sqrt(kb/kOm);
    coeff = (1i^(n + 1)) * (2*n + 1)/(n*(n + 1));
    Ben = (N*Psi_1(n,q)*Psi(n,N*q) - Psi(n,q)*Psi_1(n,N*q)) / ...
        (N*Dzeta_1(n,q)*Psi(n, N*q) - Dzeta(n, q)*Psi_1(n, N*q));
    Bmn = (N*Psi(n,q)*Psi_1(n,N*q) - Psi_1(n,q)*Psi(n,N*q)) / ...
        (N*Dzeta(n,q)*Psi_1(n,N*q) - Dzeta_1(n,q)*Psi(n,N*q));
    Ben = Ben * coeff;
    Bmn = Bmn * coeff;
    return;
end
function D = dist(Tx, Fx, rho, rx)
    global theta phi;
    syms pha(theta, phi) D(theta, phi)
    pha = sin(phi)*sin(Fx)*cos(theta-Tx)+cos(phi)*cos(Fx);
    D(theta, phi) = sqrt(rho^2 + rx^2 - 2*rho*rx*pha);
end
function Ps = Psi(n,x)
    Ps = sqrt(pi*x/2)*besselj(n,x);
end
function Dz = Dzeta(n,x)
    Dz = sqrt(pi*x/2)*besselh(n,x);
end
function val = Psi_1(n,x)
    val = sqrt(x*pi/2) * (n * besselj(n-1, x) - (n - 1)*besselj(n+1,x)) / (2*n + 1);
    val = val + sqrt(pi/(2*x))*besselh(n,x)/2;
    return;
end
function val = Dzeta_1(n,x)
    val = sqrt(pi*x/2) * (besselh(n-1,x) - (n + 1)*besselh(n,x)/x);
end
function lee = lejandr(n)
    syms f(theta) 
    global theta;
    fact1 = sqrt(1 - (cos(theta))^2);
    fact2 = -1/(pow2(n)* factorial(n));
    f(theta) = ((cos(theta))^2 - 1)^n;
    fact3 = diff(f(theta), theta, n + 1);
    lee = fact1*fact3/fact2;
end
function lee = lejandr_1(n)
    syms f(theta)
    global theta;
    lee(theta) = (theta*lejandr(n) * n - lejandr(n - 1)*(n + 1))/(theta*theta - 1);
end
dsad

Best Answer

To calculate a single integral of a symbolic expression, use int().
To calculate a double integral of a symbolic expression, use nested int() calls, like int(int(f, x, a, b), y, c, d)
Related Question