Numerical integral over the whole spherical volume of a crystallite

I am new to matlab. I want to solve the above equation to get the total carrier concentration in the spherical grain volume;

below are the values of the parameters:

Total Volume of the sphere= 4/3*pi*r^3

r= 40*10^-9 %meters

Ld= 18*10^-9 %meters

phi= 20

Nd= 1e23

Answer 1

[Torsten]

Nd = 1e23;
phi = 20;
Ld = 18e-9;
R = 40e-9;
fun = @(r) 4*pi*r.^2.*(1-exp(phi + 1/6*(r/Ld).^2));
value_integral = Nd*integral(fun,0,R)

[Anil Kumar]

Thank You for the answer sir, Thank you and respect.

Now the deal is that I need to solve this integral as a part of a non linear equation as:

Where the LHS= RHS; I need to calculate the value of unknown variable phi0

LHS is the above equation containing integral;

RHS = 4*pi*R^2*Nt/(1+2*exp(Et_Ef + K*T *(phi0 + (1/6)*(r(i)/Ld)^2)^2))

I have tried the following code:

with the following parameters:

Nd = 1e24;

e = 12*1e-12;

r = 10*10^-9

R = 10*1e-9;

Et_Ef = -0.81*1.6*1e-19;

K = 1.3807e-23;

T = 600;

Nt = 1e11;

q = 1.6e-19;

Ld = sqrt(e*K*T/(q^2*Nd));

%integrand = @(r)4*pi*(Nd-(Nd.*exp(-1.*(phi0+(1/6).*(r/Ld).^2)))).*r.^2;

fun = @(r)Nd*4*pi*r.^2.*(1-exp(phi0 + 1/6*(r/Ld).^2));

LHS = Nd*integral(fun,0,R)

%LHS = integral(integrand, 0,10*10^-9 );

RHS = 4*pi*R^2*Nt/(1+2*exp(Et_Ef + K*T *(phi0 + (1/6)*(r(i)/Ld)^2)^2));

eq = LHS - RHS==0;

phi0 = solve(eq,phi0)

But I am not able to get the value of phi0 which satisfies the LHS=RHS condition?

My question is that whether the use of solve function is ok? or shall I use vpasolve fucntion?

for vpssolve fucntion, I tried but confused how to write the LHS-RHS=0 while using vpasolve?

please demonstrate with modification in this code. I shall be highly thankful to your generousity. Thank you and respect.

[Torsten]

What is r(i) ?

Why do you multiply the LHS twice with Nd ?

Use "fzero" to solve for phi0.

[Anil Kumar]

I have tried this code;

funct= @(R,phi0) (Nd*4*pi*R.^2*(1-exp(phi0 + 1/6*(r/Ld).^2)))

equation = @(phi0)integral(@(r)funct(R,phi0),0,R,'ArrayValued',1)- 4*pi*R^2*Nt/(1+2*exp(Et_Ef + K*T *(phi0 + (1/6)*(r/Ld)^2)^2));

phi0val = fzero(equation,100)

is it the optimized approach to get phi0 in this equality condition?

[Torsten]

function main
  phi0 = ...;  % Set initial guess for phi
  phisol = fzero(@fun,phi0)  
end

function res = fun(phi)
  Nd = 1e24;
  e  = 12*1e-12;
  r  = 10*10^-9
  R  = 10*1e-9;
  Et_Ef = -0.81*1.6*1e-19;
  K  = 1.3807e-23;
  T  = 600;
  Nt = 1e11;
  q  = 1.6e-19;
  Ld = sqrt(e*K*T/(q^2*Nd));
  fun = @(r)Nd*4*pi*r.^2.*(1-exp(phi + 1/6*(r/Ld).^2));
  LHS = integral(fun,0,R);
  RHS = 4*pi*R^2*Nt/(1+2*exp(Et_Ef + K*T *(phi + (1/6)*(r/Ld)^2)^2)); %what is r here ? Should it be R ?
  res = LHS - RHS;
end

[Anil Kumar]

Thank you sir! finally fzero worked for me! thank you again for your time and humbleness! following is my code:

Nd = 1e24;

N=100

e = 12*1e-12;

r = 10*10^-9

R = 10*1e-9;

Et_Ef = -0.81*1.6*1e-19;

K = 1.3807e-23;

T = 600;

Nt = linspace(1e14,1e16,N);

q = 1.6e-19;

for i= N

Ld = sqrt(e*K*T/(q^2*Nd));

funct= @(R,phi0) (4.*pi.*(Nd-Nd.*exp(-phi0)).*R.^2)

equation = @(phi0)integral(@(R)funct(R,phi0),0,R,'ArrayValued',1)-(4*pi.*R.^2.*Nt(i)./(1+2*exp((Et_Ef + K*T .*(phi0 + (1/6).*(R/Ld).^2))/(K.*T))));

phi0val = fzero(equation,10)

end

now I got different values of phi0 at different Nt; I want to plot Nt v/s phi0? I used plot (Nt, phi0) for this, it is not working. I dont know how to plot a range of output against a range of variable? Please guide me through this.This too will be a great help, Sir!

[Torsten]

  Nd = 1e24;
  N = 100
  e  = 12*1e-12;
  R  = 10*1e-9;
  Et_Ef = -0.81*1.6*1e-19;
  K  = 1.3807e-23;
  T  = 600;
  Nt = linspace(1e14,1e16,N);
  q  = 1.6e-19;
  Ld = sqrt(e*K*T/(q^2*Nd));
  for i= 1:N 
    funct = @(r,phi0) (4*pi*(Nd-Nd*exp(-phi0))*r.^2);
    equation = @(phi0)integral(@(r)funct(r,phi0),0,R)-(4*pi*R^2*Nt(i)/(1+2*exp((Et_Ef + K*T *(phi0 + (1/6)*(R/Ld)^2))/(K*T))));
    phi0val(i) = fzero(equation,10);
  end
  plot(Nt,phi0val)

[Anil Kumar]

Thank you sir! Just one more query (please bear me!!) if I have multple values of R and I want to plot(Nt, phi0val) at different R in the same plot! how do we do that sir? Thank you in advance

below is the modified code

Nd = 1e24;

N = 100

e = 12*1e-12;

r = 10*10^-9

R = [10e-9,20e-9,30e-9,40e-9];

Et_Ef = -0.81*1.6*1e-19;

K = 1.3807e-23;

T = 600;

Nt = linspace(1e14,1e16,N);

q = 1.6e-19;

Ld = sqrt(e*K*T/(q^2*Nd));

for j = 1:length(R)

for i= 1:N

funct = @(r,phi0) (4*pi*(Nd-Nd*exp(-phi0))*r.^2);

equation = @(phi0)integral(@(r)funct(r,phi0),0,R(j))-(4*pi*R(j).^2*Nt(i)/(1+2*exp((Et_Ef + K*T *(phi0 + (1/6)*(R(j)/Ld)^2))/(K*T))));

phi0val(i) = fzero(equation,10);

end

end

plot(Nt,phi0val)

[Torsten]

  Nd = 1e24;
  N = 100;
  e  = 12*1e-12;
  R  = [10e-9,20e-9,30e-9,40e-9];
  Et_Ef = -0.81*1.6*1e-19;
  K  = 1.3807e-23;
  T  = 600;
  Nt = linspace(1e14,1e16,N);
  q  = 1.6e-19;
  Ld = sqrt(e*K*T/(q^2*Nd));
  for j = 1:numel(R)  
    for i= 1:N 
      funct = @(r,phi0) (4*pi*(Nd-Nd*exp(-phi0))*r.^2);
      equation = @(phi0)integral(@(r)funct(r,phi0),0,R(j))-(4*pi*R(j)^2*Nt(i)/(1+2*exp((Et_Ef + K*T *(phi0 + (1/6)*(R(j)/Ld)^2))/(K*T))));
      phi0val(j,i) = fzero(equation,10);
    end
  end
  plot(Nt,phi0val(1,:),Nt,phi0val(2,:),Nt,phi0val(3,:),Nt,phi0val(4,:))

[Anil Kumar]

Thanks Sir! I learnt a lot through these queries! Thanks for everything!

[Anil Kumar]

I have another issue. I want to integrate "funneff" and want to plot neff against Nt! Problem is that I am not able to integrate the function "funneff". Also how to plot neff v/s Nt ? below is my code. Thanks in advance!

Nd = 1e24;
  N = 100;
  e  = 12*1e-12;
 % R  = [10e-9,20e-9,30e-9,40e-9];
  R  = [1e-9,2e-9,3e-9,4e-9,5e-9,6e-9,7e-9,8e-9];
  Et_Ef = -1.21*1.6*1e-19;
  K  = 1.3807e-23;
  T  = 600;
  Nt = linspace(1e13,1e16,N);
  q  = 1.6e-19;
  Ld = sqrt(e*K*T/(q^2*Nd));
  for j = 1:numel(R)  
    for i= 1:N 
      funct = @(R,phi0) (4*pi*(Nd-Nd*exp(-phi0))*R.^2);
      equation = @(phi0)integral(@(R)funct(R,phi0),0,R(j))-(4*pi*R(j)^2*Nt(i)/(1+2*exp((Et_Ef + K*T *(phi0 + (1/6)*(R(j)/Ld)^2))/(K*T))));
      %equation = @(phi0)integral(@(r)funct(r,phi0),0,R)-(4*pi*R^2*Nt(i)/(1+2*exp((Et_Ef + K*T *(phi0 + (1/6)*(R/Ld)^2))/(K*T))));
      phi0val(j,i) = fzero(equation,40);
      %phi0val(i) = fzero(equation,40)
    end
  end
% plot(Nt,phi0val(1,:))
% hold on
% plot(Nt,phi0val(2,:))
% hold on
% plot(Nt,phi0val(3,:))
% hold on
% plot(Nt,phi0val(4,:))
% hold on
% plot(Nt,phi0val(5,:))
% hold on
% plot(Nt,phi0val(6,:))
% hold on
% plot(Nt,phi0val(7,:))
% hold on
% plot(Nt,phi0val(8,:))
for j = 1:numel(R) 
    for i=1:N
funneff= @(R)1./(4/3.*pi.*R(j).^3).*(Nd.*exp(-phi0val(j,i))).*R(j).^2
neff(j,i)= integral(funneff,0,R(j),'Arrayvalued',true)
    end
end
plot(neff(1,:),phi0val(1,:))
hold on
plot(neff(2,:),phi0val(2,:))
hold on
plot(neff(3,:),phi0val(3,:))
hold on
plot(neff(4,:),phi0val(4,:))
hold on
plot(neff(5,:),phi0val(5,:))
hold on
plot(neff(6,:),phi0val(6,:))
 hold on
plot(neff(7,:),phi0val(7,:))
hold on
 plot(neff(8,:),phi0val(8,:))

Answer 2

[darova]

Suggestion:

Use integral or trapz

[Anil Kumar]

Thank you sir for the suggestion.