How can i calculate double integral?

Question

Igor Arkhandeev 2020-5-17

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/526053-how-can-i-calculate-double-integral

评论： Walter Roberson 2025-6-5

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2020-5-17

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/526053-how-can-i-calculate-double-integral#answer_432906

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)

10 个评论
显示 8更早的评论隐藏 8更早的评论

Igor Arkhandeev 2020-5-17

在 MATLAB Online 中打开

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 2020-5-17

编辑：Walter Roberson 2020-5-17

dint.zip

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.

Walter Roberson 2020-5-18

在 MATLAB Online 中打开

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 2020-5-18

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!

请先登录，再进行评论。

Answer 2

Yamil Antonio 2025-6-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/526053-how-can-i-calculate-double-integral#answer_1566017

在 MATLAB Online 中打开

% Trabajo de Álgebra Lineal - Punto 2
% Autor: [Tu nombre]
% Fecha: [Fecha actual]
% Definimos la matriz identidad I y la matriz de unos J
I = eye(3);    % Matriz identidad 3x3
J = ones(3);   % Matriz de unos 3x3
% Definimos A y B según el enunciado
A = (2*I - 3*J)';    % A = transpuesta de (2I - 3J)
B = 2*I + J;         % B = (3I + J) - I = 2I + J
% Punto b: determinantes de A y B
detA = det(A);
detB = det(B);
disp('Determinante de A:');
Determinante de A:
disp(detA);
   -28
disp('Determinante de B:');
Determinante de B:
disp(detB);
    20
% Punto c: X = A transpuesta + 2B
X = A' + 2*B;
disp('Matriz X = A^T + 2B:');
Matriz X = A^T + 2B:
disp(X);
     5    -1    -1
    -1     5    -1
    -1    -1     5
% Punto d: Y = -2A + B transpuesta
Y = -2*A + B';
disp('Matriz Y = -2A + B^T:');
Matriz Y = -2A + B^T:
disp(Y);
     5     7     7
     7     5     7
     7     7     5

2 个评论
显示无隐藏无

Steven Lord 2025-6-5

This doesn't appear to have anything to do with integration, let alone computing a double integral. Did you mean to ask this as a new question? If so use the Ask link near the top of the page (and explain what difficulty you're experiencing or question you're trying to answer) to start a new thread.

Walter Roberson 2025-6-5

I do not understand how this answer fits with the original question ???

请先登录，再进行评论。

How can i calculate double integral?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论
显示 8更早的评论隐藏 8更早的评论

更多回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

How can i calculate double integral?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

10 个评论 显示 8更早的评论隐藏 8更早的评论

更多回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

10 个评论
显示 8更早的评论隐藏 8更早的评论

2 个评论
显示无隐藏无