Jump to content
  • Advertisement
Sign in to follow this  
theo2008

problem with equation in Matlab

This topic is 3663 days old which is more than the 365 day threshold we allow for new replies. Please post a new topic.

If you intended to correct an error in the post then please contact us.

Recommended Posts

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.

Share this post


Link to post
Share on other sites
Advertisement
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.

Share this post


Link to post
Share on other sites
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 :-(

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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...:-)

Share this post


Link to post
Share on other sites
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.



Share this post


Link to post
Share on other sites
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

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
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.

Share this post


Link to post
Share on other sites
Sign in to follow this  

  • Advertisement
×

Important Information

By using GameDev.net, you agree to our community Guidelines, Terms of Use, and Privacy Policy.

GameDev.net is your game development community. Create an account for your GameDev Portfolio and participate in the largest developer community in the games industry.

Sign me up!