The problem with calculating multidimensional integrals

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');

Answer 1

[Walter Roberson]

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.

[Dmitry]

Walter Roberson, That is, I have to specify it everywhere except Int_SW_inear?

[Walter Roberson]

Yes.

[Dmitry]

Walter Roberson, I did as you said and the code works. But the running time of the code is too long. Please tell me, is there any more optimal solution to this problem?

[Steven Lord]

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.

[Torsten]

@Steven Lord

https://uk.mathworks.com/matlabcentral/answers/2086623-how-can-i-calculate-my-expression-with-double-integration?s_tid=srchtitle