How can i calculate double integral?

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

Answer 1

[Walter Roberson]

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)

[Igor Arkhandeev]

Thanks for your answer, but I still can't to resolve this question. At the moment I ran code, but processing time more than 2 hour! And it still processing. I try use int(int...), ...) for inner function, but in result I get following:

int(sin(1)*(1 - cos(theta)^2)^(3/2)*(2825850565061599030300850599472384/441912275208326705232454276919945 - 13651866580585377058297446534902208i/441912275208326705232454276919945) - sin(1)*cos(theta)^2*(1 - cos(theta)^2)^(1/2)*(2825850565061599030300850599472384/441912275208326705232454276919945 - 13651866580585377058297446534902208i/441912275208326705232454276919945) - sin(1)*cos(theta)^3*sin(theta)*(1 - cos(theta)^2)^(1/2)*(7058540003372008889499610343364993024/441912275208326705232454276919945 + 1900130304961971074321645499861123072i/441912275208326705232454276919945) + sin(1)*cos(theta)*sin(theta)*(1 - cos(theta)^2)^(3/2)*(2352846667790669629833203447788331008/88382455041665341046490855383989 + 633376768320657024773881833287041024i/88382455041665341046490855383989), theta, 7029203256864545/4503599627370496, pi/2)

So, have you any idea, what's wrong? This is code for understanding. Other function attached at the top of topic.

lim2 = pi/2 - atan(a/dd);
%     innerI = @(theta, phi) abs(Escr(nmax) * gradFk()) - abs(Esct(nmax) * gradFk()) * ...
%         (-dd*cos(theta)/(sin(theta) * sin(theta)));
    innerI = @(theta, phi) Escr(nmax);
    q1 = @(theta) int(innerI, phi, 1, 2*pi);
    q2 = int(q1, theta, lim2, pi/2);
    disp(q2);

[Walter Roberson]

If it helps to have an approximation then switch to vpaintegral().

[Walter Roberson]

What are some appropriate values for f, and for nmax ?

Are you expecting the result to be complex valued?

[Walter Roberson]

I have attached revised code. The main routine is driver.m .

I think I ended up revising every function

• I got rid of the global variables
• I vectorized
• I forced use of symbolic constants to try to make more sense of the weird numbers that were showing up.

If you look at driver.m you will see two subplots. In one of the two, I use vpaintegral() of q1 to produce a symbolic expression, and then I substitute in a range of f values and plot the result. In the other of the two, I use matlabFunction to generate a numeric routine, and use integral() on it for a range of f values, and plot the results.

In each case, the result is complex-valued. It appears that the real and the imaginary part both oscillate above and below 0, so in theory there are individual real-valued points here and there.

Because my re-writes are so extensive, I recommend that you do not put this code into the same directory your code was already in; you will want your old code to refer back to.

You will see a number of references to params.Q() as a function. That function converts numeric floating point value into symbolic rationals. For reasons of numeric accuracy, converting something along the lines of x.xxxxxEyyy is better done as converting x.xxxxx to rational and multiplying it by 10^yyy.

[Igor Arkhandeev]

Thanks for your big work! I try this set of programms and it get interesting result! Unfortunately, I still can't to get result for full expression:

abs(Escr(nmax, params) * gradFk(params)) - abs(Esct(nmax, params) * gradFk(params)) * (-params.dd*cos(theta)/(sin(theta) * sin(theta)));

[Walter Roberson]

Somehow I had a typing mistake in one of the functions that escaped my testing. Corrected version is

function lee = lejandr_1(n, params)
    syms theta
    Qn = params.Q(n);
    lee(theta) = (theta*lejandr(Qn, params) .* n - lejandr(Qn - 1, params) .* (Qn + 1))./(theta*theta - 1);
end

[Igor Arkhandeev]

Yes, I see this mistake earlier. But my code can't finish.

[Igor Arkhandeev]

I find some mistakes in mathematics, which is updated version.

[Walter Roberson]

I ran the below with 250 instead of 25 points. It took 3218 seconds to compute. Looking at the results, I think you could get a very reasonable reproduction with 20 to 25 points.

syms f
nmax = 5;

params = initialize_params(f);

lim2 = params.pi/2 - atan(params.a/params.dd);

syms theta phi
innerI = abs(Escr(nmax, params) .* gradFk(params)) - abs(Esct(nmax, params) .* gradFk(params)) .* ...
        (-params.dd .* cos(theta) ./ (sin(theta) .* sin(theta)));

InnerIF = matlabFunction(innerI, 'vars', {f, theta, phi});

F = linspace(2e8,1e10,25);

tic;
Q1F = @(f,theta) integral(@(phi) InnerIF(f,theta,phi), 1, 2*pi, 'arrayvalued', true);
lim2_d = double(lim2);
Q2N = integral(@(theta) Q1F(F,theta), lim2_d, pi/2, 'arrayvalued', true);
t2 = toc;
plot(F,Q2N)
title('integral() approach');
fprintf('integral() approach took %g seconds\n', t2);

[Igor Arkhandeev]

Hi again! I reviewed the math part and found an error in my program. Thank you so much for your patience and the time given to me!