For Loop using two variables

I am new to MATLAB and I am trying to use a for loop using two variables. The question is: Generate a MATLAB program to compute and plot the Fermi function, f(E), and 1- f(E) versus ΔE = E-Ef for values of ΔE that is over the range of -0.5eV ≤ ΔE ≤ 0.5eV for varying temperature settings where Temperature = 150, 250, 350, 450 and 550K.

The code I have written is:

k = 1.3806488*10^-23; %boltzman constant
format shortEng
for dE = linspace (-0.5, 0.5, 5), T1 = 150:100:550
  fE = 1/(1+exp(dE/(k*T1))) %Fermi Function
  fE2 = 1-fE
end

However, it its only displaying T1

Answer 1

[Andrew Newell]

You need to use a couple of nested loops:

k = 1.3806488*10^-23; %boltzman constant
format shortEng
for dE = linspace (-0.5, 0.5, 5)
  for T1 = 150:100:550
      fE = 1/(1+exp(dE/(k*T1))) %Fermi Function
      fE2 = 1-fE
  end
end

For plotting these results, however, you'll need to store the values in a matrix, like this:

dE = linspace(-0.5,0.5,5);
T1 = 150:100:550;
fE = zeros(length(dE),length(T1)); % initializing fE speeds up the calculation
for i=1:length(dE)
  for j=1:length(T1)
      fE(i,j) = 1/(1+exp(dE/(k*T1))); %Fermi Function
  end
end
fE2 = 1-fE; % No need to put this in the loop

[Emile]

Thanks a lot🙏

[Angelico Manhique]

Maybe you can also simply use one FOR loop:

k = 1.3806488*10^-23; %Boltzman constant
deltaE = linspace(-0.5,0.5,5);
Temp = 150:100:550;
fE = 0;

for i = 1:length(Temp)
    fE(i) = 1 / (1+exp(deltaE(i)/(k*Temp(i)))); %Fermi Function
end

fE2 = 1-fE; %Calculating the second function 

plot(deltaE,fE,deltaE,fE2) %Plotting the results
legend('f(E)','1-f(E)')
xlabel('eV')

[Torsten]

I think fE is meant as a matrix with

fE(i,j) = 1 / (1+exp(deltaE(i)/(k*Temp(j))))

[DGM]

Consider:

% inputs
k = 1.3806488*10^-23; %boltzman constant
dE = linspace(-0.5,0.5,5);
T1 = 150:100:550;

% you could use loops
fE = zeros(length(dE),length(T1)); % initializing fE speeds up the calculation
for i = 1:length(dE)
    for j = 1:length(T1)
        fE(i,j) = 1/(1 + exp(dE(i)/(k*T1(j)))); % Fermi Function
    end
end
fE2_0 = 1-fE % No need to put this in the loop

fE2_0 = 5×5

0 0 0 0 0 0 0 0 0 0 0.5000 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% or just calculate the entire output at once
% assuming both inputs are row vectors to begin with
fE2_1 = 1 - 1./(1 + exp(dE.'./(k*T1)))

fE2_1 = 5×5

0 0 0 0 0 0 0 0 0 0 0.5000 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

% that wouldn't have been an option when this was posted.
% prior to R2016b, this is how it would be done
fE2_2 = 1 - 1./(1 + exp(bsxfun(@rdivide,dE.',k*T1)))

fE2_2 = 5×5

0 0 0 0 0 0 0 0 0 0 0.5000 0.5000 0.5000 0.5000 0.5000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Note that this looks odd because these are tiny numbers (actually several orders of magnitude smaller than floating point resolution for the numerator).

(k*T1)

ans = 1×5

1.0e-20 * 0.2071 0.3452 0.4832 0.6213 0.7594

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

As a consequence, these are huge numbers being fed to an exponential function

dE.'./(k*T1)

ans = 5×5

1.0e+20 * -2.4143 -1.4486 -1.0347 -0.8048 -0.6585 -1.2072 -0.7243 -0.5174 -0.4024 -0.3292 0 0 0 0 0 1.2072 0.7243 0.5174 0.4024 0.3292 2.4143 1.4486 1.0347 0.8048 0.6585

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

... and so we're excessively sensitive to dE and we lose all sensitivity to T1:

exp(dE.'./(k*T1))

ans = 5×5

0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 Inf Inf Inf Inf Inf Inf Inf Inf Inf Inf

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

You should probably want to rethink how you want to solve a problem like this. Naive floating point numeric calculations might not be the way to go.