how can I calculate my expression with double integration?

how can I calculate my expression with double integration?

at startup, I get an error:

Operator '.*' is not supported for operands of type 'function_handle'.

Error in sw_plt2 (line 32)

s_w_func = ( sin(k_pl.*z) .* sin(k_pl.*(z-D)) ) ./ (k_pl .* sin(k_pl*D));

My code:

z = 5:0.1:100;
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 = 10e28;
v_f = 1.57e6;
l = 1.5e-10;
tau = l / v_f;
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;
w_n = @(n) t*(2*n+1);
n_max = 49;
cuA = m / (pi*h_bar*N_0*tau);
w_z = zeros(size(z));
coef_K = cuA / (2*pi*ksi);
coef_w = (cuA * ksi * A) / (4*pi);
%syms e;
wt_sum = wt(1.2, 49);
k_pl = @(e) k0*sqrt( e + wt_sum); % k_+ off dim
s_w_func = ( sin(k_pl.*z) .* sin(k_pl.*(z-D)) ) ./ (k_pl .* sin(k_pl*D));

Operator '.*' is not supported for operands of type 'function_handle'.

S_w = imag( @(e) integral(s_w_func), 0, E ); % off dim
K = @(e, z) k0*sqrt(e + 1i*wt_sum - coef_K * S_w);
inner_fun = 1/K .* ( cos(@(z_)integral(K, z-D, z)) ./ sin(@(z_)integral(K, 0, D)) - cot(@(z_)integral(K, 0, D)));
W_res = wt_sum - 1i*coef_w .* (@(e) integral(inner_fun, 0, E)); % off-dim, norm on pi*k_b*Tc
figure;
plot(z, W_res);
xlabel('z');
ylabel('w_res(z)');
title('График w_res(z)');
grid on;
function result = wt(t, n)
result = 0;
for i = 1:n
    result = result + t * (2 * i + 1);
end
end

[Dyuman Joshi]

s_w_func = ( sin(k_pl.*z) .* sin(k_pl.*(z-D)) ) ./ (k_pl .* sin(k_pl*D));

What is this line of code supposed to do?

Please share the expression in mathematical form so we can provide appropriate suggestions.

[Dmitry]

Dyuman Joshi,

I need plot(z, W_res)

As for "s_w_func", I want to integrate this by "e".

This is the expression I want to calculate:

Answer 1

[Torsten]

Maybe you mean this, but I'm not sure:

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 = 10e28;
v_f = 1.57e6;
l = 1.5e-10;
tau = l / v_f;
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;
w_n = @(n) t*(2*n+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 = wt(1.2, 49);
k_pl = @(e) k0*sqrt( e + wt_sum); % k_+ off dim
s_w_func = @(e,z)( sin(k_pl(e).*z) .* sin(k_pl(e).*(z-D)) ) ./ (k_pl(e) .* sin(k_pl(e)*D));
S_w = @(z) imag(integral(@(e)s_w_func(e,z), 0, E,'ArrayValued',1 )); % off dim
K = @(e, z) k0*sqrt(e + 1i*wt_sum - coef_K * S_w(z));
inner_fun = @(e,z) 1/K(e,z) .* ( cos(integral(@(z_)K(e,z_), z-D, z)) ./ sin(integral(@(z_)K(e,z_), 0, D)) - cot(integral(@(z_)K(e,z_), 0, D)));
W_res = @(z)wt_sum - 1i*coef_w .* integral(@(e)inner_fun(e,z), 0, E,'ArrayValued',1); % off-dim, norm on pi*k_b*Tc

function result = wt(t, n)
result = 0;
for i = 1:n
    result = result + t * (2 * i + 1);
end
end

[Dmitry]

Torsten,

I need plot(z, W_res)

This is the expression I want to calculate:

[Torsten]

Evaluate W_res with the numerical values of your z-vector (one by one) and plot.

But integration took veeeeeeeeeeery long already for one single value for z.

[Dmitry]

Torsten,

How can I evaluate the value of W_res?

[Torsten]

First try for a single value:

W_res(1)

If this works, try for your array:

z = 5:0.1:100;
W_resnum = arrayfun(@(z) W_res(z),z)

But my guess is that the denominators in your functions will have zeros in the region of integration.

[Dmitry]

How can I make these calculations faster?

Answer 2

[John D'Errico]

The important takeaway to use is, you CANNOT perform arithmetic operations between a pair of function handles. Maybe you think you should be able to, but you CANNOT. Yours is a common error we see by new users. Language syntax is irrefutable and inescapable.

HOWEVER, you can perform those same operations on the RESULT of a function handle. For example:

f1 = @(x) x + 1;
f2 = @(x) x.^2;

f3 = @(x) f1(x).*f2(x);
f3(5)

ans = 150

By way of comparison, see what happens when you try to multiply those function handles (I must do this at the end, as MATLAB in Answers will not allow me to do anything else after an error occurs.)

f1*f2

Operator '*' is not supported for operands of type 'function_handle'.

Do you see this is the same error you got? The solution is as I showed above. You need to operate on the RESULT of the function handles, not on the function handles themselves.