How can i calculate double integral?

14 ビュー (過去 30 日間)
Igor Arkhandeev
Igor Arkhandeev 2020 年 5 月 17 日
コメント済み: Igor Arkhandeev 2020 年 5 月 18 日
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

採用された回答

Walter Roberson
Walter Roberson 2020 年 5 月 17 日
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 件のコメント
Igor Arkhandeev
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!

サインインしてコメントする。

その他の回答 (0 件)

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by