How can i calculate double integral?
3 views (last 30 days)
Show older comments
Igor Arkhandeev
on 17 May 2020
Commented: Walter Roberson
on 5 Jun 2025
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 Comments
Accepted Answer
Walter Roberson
on 17 May 2020
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 Comments
Walter Roberson
on 18 May 2020
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);
More Answers (1)
Yamil Antonio
on 5 Jun 2025
% 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:');
disp(detA);
disp('Determinante de B:');
disp(detB);
% Punto c: X = A transpuesta + 2B
X = A' + 2*B;
disp('Matriz X = A^T + 2B:');
disp(X);
% Punto d: Y = -2A + B transpuesta
Y = -2*A + B';
disp('Matriz Y = -2A + B^T:');
disp(Y);
2 Comments
Steven Lord
on 5 Jun 2025
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
on 5 Jun 2025
I do not understand how this answer fits with the original question ???
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!