The problem with calculating multidimensional integrals

2 次查看（过去 30 天）
Dmitry
Dmitry 2024-2-28
评论： Torsten 2024-2-29
Ran in:
The problem with calculating multidimensional integrals.
I'm trying to calculate my integrals and plot a graph, but it gives me an error when I start: "Arrays have incompatible sizes for this operation."
How can I solve this problem so that the code works?
My code:
z = 5:0.1:70;
D = 100;
e_f = 1.88e-18;
m = 9.11e-31;
k_b = 1.38e-23;
h_bar = 1.05e-34;
ksi = 1e-9;
t = 1;
A = 2 * m / h_bar^2;
k_f = sqrt(e_f / (pi * k_b * 1.2));
N_0 = 1e28;
v_f = 1.57e6;
l = 1.5e-10;
tau = h_bar / (3.056 * 2 * pi * k_b * 1.2);
gamma = h_bar / (2 * pi * k_b * 1.2 * tau);
E = e_f / (pi * k_b * 1.2);
norm = 2 * m / h_bar^2;
k0 = (ksi / h_bar) * sqrt(2 .* m .* pi .* k_b .* 1.2);
coef = 1;
n_max = 49;
cuA = m / (pi * h_bar * N_0 * tau);
coef_K = cuA / (2 * pi * ksi);
coef_w = (cuA * ksi * A) / (4 * pi);
wt_sum = 2499;
wz = zeros(size(z));
cos_K_E_res = zeros(size(z));
Wn_2 = zeros(size(z));
for i = 1:length(z)
k_plus = @(e) k0*sqrt( e + 1i*wt_sum );
Sw_inear = @(e, zz_) (sin(k_plus(e) .* zz_) .* sin(k_plus(e) .* (zz_ - D))) ./ (k_plus(e) .* sin(k_plus(e) *D));
Int_SW_inear = @(zz_) integral(@(e) Sw_inear(e, zz_), 0, E);
K_z = @(e, zz_) k0*sqrt( e + 1i*wt_sum - coef_K * imag(Int_SW_inear(zz_)));
trigl_funcs = @(e) cos( integral(@(zz_) K_z(e, zz_), D-z(i) , z(i)) ) ./ sin( integral(@(zz_) K_z(e, zz_), 0, D) ) - cot( integral(@(zz_) K_z(e, zz_), 0, D) );
w_inear = @(e) trigl_funcs(e) ./ K_z(e, z(i));
%Wn_2(i) = wt_sum + 1i*coef_w * integral(w_inear, 0, E, 'ArrayValued', 1);
%Wn_2(i) = wt_sum + 1i * coef_w * integral(@(e) w_inear, 0, E,'ArrayValued', 1);% %%%%%%%%%
Wn_2(i) = wt_sum + 1i * coef_w * integral(@(e) w_inear(e), 0, E, 'ArrayValued', 1);
%Inear_cos_K = integral(@(zz) K(E, zz), z(i) - D, z(i));
%cos_K_E_res(i) = integral(@(e) Inear_cos_K, 0, E, 'ArrayValued', 1);
end
Arrays have incompatible sizes for this operation.

Error in solution>@(e,zz_)(sin(k_plus(e).*zz_).*sin(k_plus(e).*(zz_-D)))./(k_plus(e).*sin(k_plus(e)*D)) (line 30)
Sw_inear = @(e, zz_) (sin(k_plus(e) .* zz_) .* sin(k_plus(e) .* (zz_ - D))) ./ (k_plus(e) .* sin(k_plus(e) *D));

Error in solution>@(e)Sw_inear(e,zz_) (line 31)
Int_SW_inear = @(zz_) integral(@(e) Sw_inear(e, zz_), 0, E);

Error in integralCalc/iterateScalarValued (line 323)
fx = FUN(t).*w;

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in solution>@(zz_)integral(@(e)Sw_inear(e,zz_),0,E) (line 31)
Int_SW_inear = @(zz_) integral(@(e) Sw_inear(e, zz_), 0, E);

Error in solution>@(e,zz_)k0*sqrt(e+1i*wt_sum-coef_K*imag(Int_SW_inear(zz_))) (line 32)
K_z = @(e, zz_) k0*sqrt( e + 1i*wt_sum - coef_K * imag(Int_SW_inear(zz_)));

Error in solution>@(zz_)K_z(e,zz_) (line 34)
trigl_funcs = @(e) cos( integral(@(zz_) K_z(e, zz_), D-z(i) , z(i)) ) ./ sin( integral(@(zz_) K_z(e, zz_), 0, D) ) - cot( integral(@(zz_) K_z(e, zz_), 0, D) );

Error in integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

Error in integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);

Error in solution>@(e)cos(integral(@(zz_)K_z(e,zz_),D-z(i),z(i)))./sin(integral(@(zz_)K_z(e,zz_),0,D))-cot(integral(@(zz_)K_z(e,zz_),0,D)) (line 34)
trigl_funcs = @(e) cos( integral(@(zz_) K_z(e, zz_), D-z(i) , z(i)) ) ./ sin( integral(@(zz_) K_z(e, zz_), 0, D) ) - cot( integral(@(zz_) K_z(e, zz_), 0, D) );

Error in solution>@(e)trigl_funcs(e)./K_z(e,z(i)) (line 36)
w_inear = @(e) trigl_funcs(e) ./ K_z(e, z(i));

Error in solution>@(e)w_inear(e) (line 39)
Wn_2(i) = wt_sum + 1i * coef_w * integral(@(e) w_inear(e), 0, E, 'ArrayValued', 1);

Error in integralCalc/iterateArrayValued (line 156)
fxj = FUN(t(1)).*w(1);

Error in integralCalc/vadapt (line 130)
[q,errbnd] = iterateArrayValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral (line 87)
Q = integralCalc(fun,a,b,opstruct);
plot(z, real(Wn_2), '-r', z, imag(Wn_2), '-b');
legend('real part', 'imag part');

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

回答（1 个）

Walter Roberson
Walter Roberson 2024-2-28
Each layer of integral() that you call passes in a variable-length vector of values.
trigl_funcs = @(e) cos( integral(@(zz_) K_z(e, zz_), D-z(i) , z(i)) ) ./ sin( integral(@(zz_) K_z(e, zz_), 0, D) ) - cot( integral(@(zz_) K_z(e, zz_), 0, D) );
A variable length vector of values will get passed as zz_ to K_z
K_z = @(e, zz_) k0*sqrt( e + 1i*wt_sum - coef_K * imag(Int_SW_inear(zz_)));
which will pass the variable-length vector to Int_SW_inear
Int_SW_inear = @(zz_) integral(@(e) Sw_inear(e, zz_), 0, E);
which will pass it to Sw_inear, along with a different variable-length vector of values, e
Sw_inear = @(e, zz_) (sin(k_plus(e) .* zz_) .* sin(k_plus(e) .* (zz_ - D))) ./ (k_plus(e) .* sin(k_plus(e) *D));
which will try to .* together k_plus(e) and zz_ -- but e and zz_ are different lengths vectors of values.
To get around this problem, you need to specify 'ArrayValued', 1 on every layer of integral other than the first.
  5 个评论
显示 3更早的评论
Steven Lord
Steven Lord 2024-2-29
I ran your code (with ArrayValued set to true in all the integral calls) in the Profiler for just the first two values of z. It calls integral 136,202 times taking just over 8 seconds on my machine. The most expensive operation is one of the helpers for integral, but most of the time required by that helper is evaluating your @(e) Sw_inear(e, zz_) function 20,295,000 times (taking a total time of almost 3 minutes.) But that's not the function that gets called most frequently; that would be @(e) k0*sqrt(e+1i*wt_sum) which is run 81,180,000 times!
And when I look at the figure that gets created and at the computed values, I'm a little skeptical of the values being computed.
>> Wn_2(1:2)
ans =
Column 1
-3.67384930932042e+19 + 6.82661560825731e+18i
Column 2
-3.67391096858457e+19 + 6.82747977875181e+18i
You've got a lot of code there so it's difficult to understand what exactly you're trying to compute. I suspect the fact that your constants span many orders of magnitude, from h_bar around 1e-34 to N_0 at 1e28, also contributes to the values of Wn_2 having components on the order of 1e18 or 1e19.
Perhaps if you showed us the mathematical form of the equation you're trying to implement to calculate Wn_2 we may be able to offer some suggestions for optimizations. You could use the LaTeX button in the MATLAB Answers Editor (the Sigma sign in the Insert section of the toolstrip) to enter it using LaTeX.
At a glance, though, one (potentially small) optimization might be to define trigl_funcs as a local function rather than an anonymous function. Then you could compute the integral that you call sin and cot on once rather than computing it twice, once to take the sine and once to take the cotangent.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by