Plotting the function by the points that need to be determined

It is necessary to plot F(t) by points. The function F(t) is a sum from 0 to nD(t)(nD is the upper limit of the sum) which depends on 't', i.e. I have an array 't' and an array nD(t) is formed from it, it contains 20 values [496, 248, 165, ...], the first point will be the final sum of F(t) with an upper limit of 496, the first point will be the final sum of F(t) with an upper limit of 248, etc., it is necessary to plot F(t) at these 20 points.

My code:

%% initial conditions

global d k0 h_bar ksi m E;

Ef = 2.77*10^3;

Kb = physconst('boltzmann'); % 1.38*10^(-23)

T = 0.12:0.24:6.4;

m = 9.1093837*10^(-31);

Tc = 1.2;

%t = T./Tc;

t = 0.1:0.1:2;

nD = floor(375./(2.*pi.*t.*1.2) - 0.5);

D = 10^(-8); % толщина пленки

ksi = 10^(-9);

%d = D/ksi;

d = 1000;

E = Ef/(pi*Kb*Tc);

h_bar = (1.0545726*10^(-34));

k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);

C_2 = 0;

for n = 0:49

C_2 = C_2 + (1/(2.*n+1)).*k0.*real(sqrt(3601+1i.*(2.*n+1))-((1+1i)./sqrt(2)).*sqrt(2.*n+1)); % константа

end

%% calculation

F = f_calc(t,nD);

plot(t,F, '-r');

%% F(t)

function F = f_calc(t,nD)

global d k0 h_bar ksi m;

F = 0;

for i = 1:20

n = nD(1,i);

F = F + 1/(2*n+1).*(k0.*real(((f_p1(n,t)-f_p2(n,t))./2))+(f_arg_2(n,t)-f_arg_1(n,t))./d);

end

F = -F;

%F = -(1/d).*F;

%F = F - C_2;

end

function p1 = f_p1(n,t)

p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));

end

function p2 = f_p2(n,t)

global E;

p2 = sqrt(3601+1i.*t.*(2.*n+1));

end

function n_lg = f_lg(n,t)

global d k0;

arg_of_lg = (1+exp(-1i*d*k0.*f_p1(n,t)))/(1+exp(-1i*d*k0.*f_p2(n,t)));

n_lg = log(abs(arg_of_lg));

end

function arg_1 = f_arg_1(n,t)

global d k0;

arg_1 = angle(1+exp(-1i*d*k0.*f_p1(n,t)));

end

function arg_2 = f_arg_2(n,t)

global d k0;

arg_2 = angle(1+exp(-1i*d*k0.*f_p2(n,t)));

end

[Torsten]

So you expect F to be a matrix of size numel(t) x numel(nD) where element (i,j) is the sum for t(i), taken from 0 to nD(j) ?

[Dmitry]

Yes, you are right.

Answer 1

[Torsten]

Maybe something like this ?

%% initial conditions
global d k0 h_bar ksi m E;
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
T = 0.12:0.24:6.4;
m = 9.1093837*10^(-31);
Tc = 1.2;
%t = T./Tc;
t = 0.1:0.1:2;
nD = floor(375./(2.*pi.*t.*1.2) - 0.5);
D = 10^(-8); % толщина пленки
ksi = 10^(-9);
%d = D/ksi;
d = 1000;
E = Ef/(pi*Kb*Tc);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 0;
for n = 0:49
    C_2 = C_2 + (1/(2.*n+1)).*k0.*real(sqrt(3601+1i.*(2.*n+1))-((1+1i)./sqrt(2)).*sqrt(2.*n+1)); % константа
end
%% calculation
F = f_calc(t,nD);
hold on
plot(t,F(:,1),"Color","red");
plot(t,F(:,numel(nD)),"Color","blue");
hold off
grid on

function F = f_calc(t,nD) 
global d k0 h_bar ksi m;
F = zeros(numel(t),numel(nD));
for i = 1:numel(t)
   for j = 1:numel(nD)
      n = nD(j);  
      for k = 0:n
        F(i,j)  = F(i,j) + 1/(2*k+1).*(k0.*real(((f_p1(k,t(i))-f_p2(k,t(i)))./2))+(f_arg_2(k,t(i))-f_arg_1(k,t(i)))./d); 
      end
   end
end
F = -F;
%F = -(1/d).*F;
%F = F - C_2;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
global E;
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function n_lg = f_lg(n,t)
global d k0;
arg_of_lg = (1+exp(-1i*d*k0.*f_p1(n,t)))/(1+exp(-1i*d*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end
function arg_1 = f_arg_1(n,t)
global d k0;
arg_1 = angle(1+exp(-1i*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t)
global d k0;
arg_2 = angle(1+exp(-1i*d*k0.*f_p2(n,t)));
end

[Dmitry]

Please explain the line plot(t.',[F(:,1),F(:,numel(nD))], '-r');

[Torsten]

The command plots F(t) against t where in the calculation for F, nD(1) (1st curve) and nD(20) (2nd curve) terms have been used.

[Dmitry]

That is, there are 2 graphs on the screen, the first F(t) is upper and the second nD(t) is lower?

[Dmitry]

Why do we have 2 curves?

[Torsten]

Yes, you plot F(:,1) against t and F(:,numel(nD)) against t. This gives two curves.

Choose different colors if you want to see which curve is F(:,1) and which is F(:,numel(nD)).

I changed the graphics above - maybe you like it more like this.

[Dmitry]

I think we misunderstood each other:

We have:

Each of these sums is finite, there are 20 such sums in total, since there are twenty values in the array t[]. Therefore, our final graph will consist of these 20 points.

[Torsten]

And what's the sense that you take a different number of elements for the sum depending on t ?

Do you know in advance how fast the infinite series converges depending on t ?

See the answer below.

Answer 2

[Torsten]

%% initial conditions
global d k0 h_bar ksi m E;
Ef = 2.77*10^3; 
Kb = physconst('boltzmann'); % 1.38*10^(-23)
T = 0.12:0.24:6.4;
m = 9.1093837*10^(-31);
Tc = 1.2;
%t = T./Tc;
t = 0.1:0.1:2;
nD = floor(375./(2.*pi.*t.*1.2) - 0.5);
D = 10^(-8); % толщина пленки
ksi = 10^(-9);
%d = D/ksi;
d = 1000;
E = Ef/(pi*Kb*Tc);
h_bar = (1.0545726*10^(-34));
k0 = (ksi/h_bar)*sqrt(2.*m.*pi.*Kb.*Tc);
C_2 = 0;
for n = 0:49
    C_2 = C_2 + (1/(2.*n+1)).*k0.*real(sqrt(3601+1i.*(2.*n+1))-((1+1i)./sqrt(2)).*sqrt(2.*n+1)); % константа
end
%% calculation
F = f_calc(t,nD);
plot(t,F)
grid on

function F = f_calc(t,nD) 
global d k0 h_bar ksi m;
F = zeros(1,numel(t));
for i = 1:numel(t) 
  for k = 0:nD(i)
    F(i)  = F(i) + 1/(2*k+1).*(k0.*real(((f_p1(k,t(i))-f_p2(k,t(i)))./2))+(f_arg_2(k,t(i))-f_arg_1(k,t(i)))./d); 
  end
end
F = -F;
%F = -(1/d).*F;
%F = F - C_2;
end
function p1 = f_p1(n,t)
p1 = ((1+1i)./sqrt(2)).*sqrt(t.*(2.*n+1));
end
function p2 = f_p2(n,t)
global E;
p2 = sqrt(3601+1i.*t.*(2.*n+1));
end 
function n_lg = f_lg(n,t)
global d k0;
arg_of_lg = (1+exp(-1i*d*k0.*f_p1(n,t)))/(1+exp(-1i*d*k0.*f_p2(n,t)));
n_lg = log(abs(arg_of_lg));
end
function arg_1 = f_arg_1(n,t)
global d k0;
arg_1 = angle(1+exp(-1i*d*k0.*f_p1(n,t)));
end
function arg_2 = f_arg_2(n,t)
global d k0;
arg_2 = angle(1+exp(-1i*d*k0.*f_p2(n,t)));
end