problem with equation in Matlab
Author
Hi everyone.A have a problem with my equation and i need help.
I have this big equation:
Cdot=qAL∂/∂v∫(1/(1+(exp(E-qVdot)/kT))*(gNdot/sqrt(piM/2))
*exp(-2((E+Edot+qVdot)/M)^2)dE.
Now i want to insert this equation in Matlab(sorry,my English is not very good)
I have these values
A=12.56*10^-4;,L=0.44;,G=2;,T=300;,q=1.6*10^-19;,k=8.62*10^-5;,
Edot=0.01;,Ndot=10^8;,M=0.001;
Edot must be between 0.01-1
Ndot between 10^8-10^15
M between 0.001-0.5
The integral i want to be (0,inf)
V must be:V=-12:0.05:2;
Also,Vdot=V/L;
First of all,I want this big function Cdot to give me a result for each value of V every time.After i must create a plot(V,Cdot).
The code that i am thinking is:
First i create a function in a M-file with name integral3:
function total=integral3(vdot)
A=12.56*10^-4;
L=0.44;
G=2;
T=300;
q=1.6*10^-19;
k=8.62*10^-5;
Edot=0.01;
Ndot=10^8;
M=0.001;
de=0.01;
E=0:0.01:1000;
expe=((E+Edot+(q*vdot)/L)/M).^2;
den=1+exp((E-(q*vdot)/L)/(k*T));
static=(G*Ndot)/(sqrt(pi/2)*M);
sinar2=exp(-2*expe);
sinar1=1./den;
sinar3=sinar1.*static;
final=sinar3.*sinar2;
subtotal=sum(final)*de;
total=q*A*L*subtotal;
And after in other M-file:
cdot1(1)=integral3(-12);
cdot2(281)=integral3(2.05);
for i=2:281
vdot(i)=-12+(i-1)*0.05;
cdot1(i)=integral3(vdot(i));
cdot2(i-1)=cdot1(i);
cdot(i)=1e11*(cdot2(i)-cdot1(i))/0.05;
end
plot(vdot,cdot,'.:')
My problem is that with this loop Cdot has always the same value...What going wrong?
What can i change to have for different value of Vdot,different Cdot?
Thank you.
I think the reason that in your 'integral3' function returns the same value all the time is that the value you pass in, 'Vdot', is multiplied by 'q', which is a tiny number. So it has a very small effect on the result.
You should consider non-dimensionalizing your equation. It will help you to figure out if your equation can be simplified, and it can make your equation more solver-friendly.
Also, Matlab is pretty good at integration so maybe you should get it to solve your differential equation.
You should consider non-dimensionalizing your equation. It will help you to figure out if your equation can be simplified, and it can make your equation more solver-friendly.
Also, Matlab is pretty good at integration so maybe you should get it to solve your differential equation.
Author
thank you.But i dont know if wikipedia can help me with this problem.A have this problem 2 months now i dont know what can i change in my code to have the results i want :-(
Quote:Original post by theo2008
thank you.But i dont know if wikipedia can help me with this problem.A have this problem 2 months now i dont know what can i change in my code to have the results i want :-(
Ok, I may be able to help but I need to know some more information about the original problem. It looks like some kind of statistical mechanics flux or partition problem, but I don't really know what it is.
You have an equation for Cdot. In order to calculate that, you need to solve an integral. I need to know what the variables in that integral stand for physically.
The second exp term in the integral looks wrong. If Edot is the time derivative of E, it must have different dimensions. So the expression ((E + Edot + qVdot)/M) must have conflicting dimensions, ie they cannot be added together.
Expressions inside exp() must not have physical dimensions. Look at the first exp term, (E + qVdot)/kT. Assuming that E is energy, qVdot is some kind of linear kinetic energy, k is Boltzmann's constant, and T is temperature, you have a numerator and denominator with the same physical units. So the physical units cancel out.
Check that exp() term and then we can see if we can make any more progress on the numerics.
Author
thank you for your notes,tomorrow i'll see again my code and i'll tell you what's happen exactly...
Author
Because i made some mistakes with my equation,the correct equation is this:
i give it to you with the parameters a,b,c because i want to understand completely.
a=(1/(1+(exp(E-qVdot)/kT))
b=(gNdot/sqrt(piM/2))
c=exp(-2*((E+Edot+qVdot)/M)^2)
So the final equation is
Cdot=qAL∂/∂v∫(a*b*c)dE
Now,
a=diode area
L=level arm coefficient
q=charge
k=Boltzmann constant
g=spin degeneracy
Edot=quantum dots energy
Yesterday i make some changes in my program:
i wrote this line
cdot(i)=1e11*(cdot2(i)-cdot1(i))/0.05;
without 1e11
cdot(i)=(cdot2(i)-cdot1(i))/0.05;
Finally my program which i created is:
function total=integral3(v)
A=12.56*10^-4;
L=0.44;
G=2;
T=300;
q=1.6*10^-19;
k=8.62*10^-5;
Edot=0.01;
Ndot=10^8;
M=0.001;
de=0.01;
vdot=v/L;
E=0:0.01:1000;
expe=((E+Edot+(q*vdot))/M).^2;
den=1+exp((E-(q*vdot))/(k*T));
static=(G*Ndot)/(sqrt(pi*M/2));
sinar2=exp(-2*expe);
sinar1=1./den;
sinar3=sinar1.*static;
final=sinar3.*sinar2;
subtotal=sum(final)*de;
total=q*A*L*subtotal;
and one other M-file:
cdot1(1)=integral3(-12);
cdot2(281)=integral3(2.05);
for i=2:281
v(i)=-12+(i-1)*0.05;
cdot1(i)=integral3(v(i));
cdot2(i-1)=cdot1(i);
cdot(i)=(cdot2(i)-cdot1(i))/0.05;
end
plot(v,cdot,'.:')
I changed my code because i think that my previous code didnt be with respect to v.I 'run' again the program.Now i have a better result something like -1.9134e-015 but this result is again the same for all the values of v and again my plot is straight line.something goes wrong with my problem.But what is the question...:-)
i give it to you with the parameters a,b,c because i want to understand completely.
a=(1/(1+(exp(E-qVdot)/kT))
b=(gNdot/sqrt(piM/2))
c=exp(-2*((E+Edot+qVdot)/M)^2)
So the final equation is
Cdot=qAL∂/∂v∫(a*b*c)dE
Now,
a=diode area
L=level arm coefficient
q=charge
k=Boltzmann constant
g=spin degeneracy
Edot=quantum dots energy
Yesterday i make some changes in my program:
i wrote this line
cdot(i)=1e11*(cdot2(i)-cdot1(i))/0.05;
without 1e11
cdot(i)=(cdot2(i)-cdot1(i))/0.05;
Finally my program which i created is:
function total=integral3(v)
A=12.56*10^-4;
L=0.44;
G=2;
T=300;
q=1.6*10^-19;
k=8.62*10^-5;
Edot=0.01;
Ndot=10^8;
M=0.001;
de=0.01;
vdot=v/L;
E=0:0.01:1000;
expe=((E+Edot+(q*vdot))/M).^2;
den=1+exp((E-(q*vdot))/(k*T));
static=(G*Ndot)/(sqrt(pi*M/2));
sinar2=exp(-2*expe);
sinar1=1./den;
sinar3=sinar1.*static;
final=sinar3.*sinar2;
subtotal=sum(final)*de;
total=q*A*L*subtotal;
and one other M-file:
cdot1(1)=integral3(-12);
cdot2(281)=integral3(2.05);
for i=2:281
v(i)=-12+(i-1)*0.05;
cdot1(i)=integral3(v(i));
cdot2(i-1)=cdot1(i);
cdot(i)=(cdot2(i)-cdot1(i))/0.05;
end
plot(v,cdot,'.:')
I changed my code because i think that my previous code didnt be with respect to v.I 'run' again the program.Now i have a better result something like -1.9134e-015 but this result is again the same for all the values of v and again my plot is straight line.something goes wrong with my problem.But what is the question...:-)
Firstly, you want to make sure that you go as far with the mathematics as you can *before* you start solving the problem numerically. It will often give you a much more stable problem to solve. So, the first thing you should do is use Leibniz Rule so you don't need to calculate the derivative.
Also, what is the interval of integration? Is it really from 0 to 1000? Or is it 0 to infinity?
I'll take a look at the numerics and see if I can find anything helpful.
Also, what is the interval of integration? Is it really from 0 to 1000? Or is it 0 to infinity?
I'll take a look at the numerics and see if I can find anything helpful.
Author
my interval is 0 to inf.i used the 0-1000 because i dont know how to define the 0-inf in matlab.
hmmm,its good idea to use the leibniz rule but my derivative is with respect to v and integral with respect to E if i dont make a mistake!!I dont know if i have problem with it..
one of my friends suggest me Gauss-konrod quadrature.What do you think about this method?
http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/ref/quadgk.html&http://www.mathworks.com/cgi-bin/texis/webinator/search/?db=MSS&prox=page&rorder=750&rprox=750&rdfreq=500&rwfreq=500&rlead=250&sufs=0&order=r&is_summary_on=1&ResultCount=10&query=quadgk&submitButtonName=Search
hmmm,its good idea to use the leibniz rule but my derivative is with respect to v and integral with respect to E if i dont make a mistake!!I dont know if i have problem with it..
one of my friends suggest me Gauss-konrod quadrature.What do you think about this method?
http://www.mathworks.com/access/helpdesk/help/techdoc/index.html?/access/helpdesk/help/techdoc/ref/quadgk.html&http://www.mathworks.com/cgi-bin/texis/webinator/search/?db=MSS&prox=page&rorder=750&rprox=750&rdfreq=500&rwfreq=500&rlead=250&sufs=0&order=r&is_summary_on=1&ResultCount=10&query=quadgk&submitButtonName=Search
Quote:Original post by theo2008
my interval is 0 to inf.i used the 0-1000 because i dont know how to define the 0-inf in matlab.
hmmm,its good idea to use the leibniz rule but my derivative is with respect to v and integral with respect to E if i dont make a mistake!!I dont know if i have problem with it..
That's exactly when it should be used. Take another look at the link I posted.
Quote:Original post by theo2008
one of my friends suggest me Gauss-konrod quadrature.What do you think about this method?
That looks fine, and 'quadgk' handles infinity for you.
Essentially the issue is that qVdot is so small that it has almost no effect upon the integral.
Take a look at each of the exponential terms in the integral.
The first exponential term:
1 / ( 1 + exp( ( E - qVdot ) / kT ) )
= 1 / ( 1 + exp( -qVdot / kT ) exp( E / kT ) )
≅ 1 / ( 1 + exp( E / kT ) )
because exp( -qVdot / kT ) ≅ 1
The second exponential term:
exp( -2 ( ( E + Edot + qVdot ) / M )^2 )
≅ exp( -2 ( ( E + Edot ) / M )^2 )
because Edot is huge compared to qVdot.
If this is not what you would expect then you need to check that you have the correct equation. No numerical technique is going to be able to reliably resolve the effect of V at these scales.
Take a look at each of the exponential terms in the integral.
The first exponential term:
1 / ( 1 + exp( ( E - qVdot ) / kT ) )
= 1 / ( 1 + exp( -qVdot / kT ) exp( E / kT ) )
≅ 1 / ( 1 + exp( E / kT ) )
because exp( -qVdot / kT ) ≅ 1
The second exponential term:
exp( -2 ( ( E + Edot + qVdot ) / M )^2 )
≅ exp( -2 ( ( E + Edot ) / M )^2 )
because Edot is huge compared to qVdot.
If this is not what you would expect then you need to check that you have the correct equation. No numerical technique is going to be able to reliably resolve the effect of V at these scales.
This topic is closed to new replies.
Advertisement
Popular Topics
Advertisement
