Integrating a complex function

8 次查看(过去 30 天)
Dmitry
Dmitry 2023-10-30
评论: Star Strider 2023-10-30
I want to count such an expression:
This function depends on "D", and "k" is why I'm doing the integration. The sum is here because my complex integration limits "k1" and "k2" depend on "en". I calculate all this, but the resulting array turns out to be empty, with both the real and imaginary parts.
Thank you very much in advance!
my code:
z = 5;
D = 5:0.1: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;
k_f = (1/h_bar)*sqrt(2*m*e_f); % k_F
N_0 = k_f^3/(3*pi^2); % Here k = k_F
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);
nD = floor(375/(2*pi.*t*1.2) - 0.5);
% Pre-allocation of the array for results
result = zeros(size(D));
for j = 1:length(D)
S12 = 0;
for n = 0:nD
k1 = (1 + 1i)/sqrt(2) * sqrt(t * (2 * n + 1) + gamma); % Set the value to k1
k2 = (1 + 1i)/sqrt(2) * sqrt(E + 1i*(t * (2 * n + 1)) + 1i*gamma); % Set the value to k2
S12 = S12 + integral(@(k) (sin(k.*z) .* sin(k.*(z-D(j))) ) ./ sin(k.*D(j)), k1, k2);
end
result(j) = S12;
end
plot(D, result);
  3 个评论
Dmitry
Dmitry 2023-10-30
编辑:Dmitry 2023-10-30
Star Strider, that is, the problem is within the integration, and the rest of the code is correct?
Star Strider
Star Strider 2023-10-30
@Dmitry — I cannot determine that. I do not understand what the code is supposed to do.
I can only say that it is not correct if the results are uniformly NaN. Something is obviously wrong, because the NaN values are likely the result of the sin functions returning only values for their arguments that are integer multiples of pi. I suspect that is not what you want.

请先登录,再进行评论。

回答(1 个)

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2023-10-30
Here are two points to consider:
(1) to specify the error tolerances and make them more crude, e.g.: 'RelTol',10,'AbsTol',0.1
(2) to compute Real and Imag parts of the function in separate corresponding to real and imaginary values of k1 and k2.
z = 5;
D = 5:0.1: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;
k_f = (1/h_bar)*sqrt(2*m*e_f); % k_F
N_0 = k_f^3/(3*pi^2); % Here k = k_F
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);
nD = floor(375/(2*pi.*t*1.2) - 0.5);
%%
% Pre-allocation of the array for results
result = zeros(size(D));
for j = 1:length(D)
S12 = 0;
SI = 0;
SR = 0;
for n = 0:nD
k1 = (1 + 1i)/sqrt(2) * sqrt(t * (2 * n + 1) + gamma); % Set the value to k1
k2 = (1 + 1i)/sqrt(2) * sqrt(E + 1i*(t * (2 * n + 1)) + 1i*gamma); % Set the value to k2
k1R=real(k1); k2R=real(k2);
k1I=imag(k1); k2I=imag(k2);
SR = SR + integral(@(k) (sin(k.*z) .* sin(k.*(z-D(j))) ) ./ sin(k.*D(j)), k1R, k2R, 'RelTol',10,'AbsTol',0.1);
SI = SI + integral(@(k) (sin(k.*z) .* sin(k.*(z-D(j))) ) ./ sin(k.*D(j)), k1I, k2I, 'RelTol',10,'AbsTol',0.1);
% S12 = S12 + integral(@(k) (sin(k.*z) .* sin(k.*(z-D(j))) ) ./ sin(k.*D(j)), k1, k2, 'RelTol',10,'AbsTol',0.1);
end
R(j) = SR;
I(j) = SI;
% result(j) = S12;
end
%plot(D, result);
nexttile
plot(D, R)
ylabel('Real Part')
grid on
nexttile
plot(D, I)
ylabel('Imag part')
xlabel('D')
grid on
% If necessary, then Real and Imag parts can be combined
  1 个评论
Torsten
Torsten 2023-10-30
编辑:Torsten 2023-10-30
I get a wrong result when I apply your way of complex integration.
Did I misunderstand your code suggestion ?
% Usual definition of complex integral (path-independent since f is
% holomorphic)
f = @(z) z.^2;
k1 = 1i;
k2 = -1 - 2*1i;
k = @(t) (1-t)*k1 + t*k2;
kdot = @(t)-k1+k2;
fun = @(t) f(k(t))*kdot(t);
I1 = integral(fun,0,1)
I1 = 3.6667 + 1.0000i
% Evaluation via antiderivative
1/3*(k2^3-k1^3)
ans = 3.6667 + 1.0000i
% Your way to integrate
k1R=real(k1); k2R=real(k2);
k1I=imag(k1); k2I=imag(k2);
IR = integral(f,k1R,k2R);
II = integral(f,k1I,k2I);
I2 = IR + 1i*II
I2 = -0.3333 - 3.3333i

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

产品

Community Treasure Hunt

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

Start Hunting!

Translated by