Integrating a complex function

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

[Star Strider]

The integrand is uniformly NaN (specifically NaN + NaNi) each for ‘k1’ and ‘k2’ for all values of ‘n’. That is the result of 0/0 or Inf/Inf so the numerator and denominator are evaluating to those values in each iteration.

I have no specific suggestions on how to solve that problem.

[Dmitry]

Star Strider, that is, the problem is within the integration, and the rest of the code is correct?

[Star Strider]

@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.

Answer 1

[Sulaymon Eshkabilov]

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

[Torsten]

@Sulaymon Eshkabilov

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