The problem with calculating multidimensional integrals
2 次查看(过去 30 天)
显示 更早的评论
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
plot(z, real(Wn_2), '-r', z, imag(Wn_2), '-b');
legend('real part', 'imag part');
0 个评论
回答(1 个)
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 个评论
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 中查找有关 Switches and Breakers 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!