How can i calculate double integral?

3 views (last 30 days)
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

Accepted Answer

Walter Roberson
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
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);
Igor Arkhandeev
Igor Arkhandeev on 18 May 2020
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!

Sign in to comment.

More Answers (1)

Yamil Antonio
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:');
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 Comments
Steven Lord
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
Walter Roberson on 5 Jun 2025
I do not understand how this answer fits with the original question ???

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!