2
$\begingroup$

I am solving a second order differential equation for n[x]. I solved it with both DSolve and NDSolve. A plot of n as a function of x looks great. Than I want to calculate a gradient of n[x] at a point x=0 and plot it as a function of V. I don't know how to execute this. I calculated gradient at every point and plot it. I get expected shape of the curve and there is also V on x axis, but the values of y axis are to high. Here is the code:

eq = Gx + D[Dif*\!\(
\*SubscriptBox[\(\[PartialD]\), \(x\)]\(n[x]\)\), x] - 
    Kr (n[x] - n0) == 0;
bcs = {n'[d] == 0, n[0] == n0*Exp[(e \[Alpha] V)/(kb T)]};
\[Alpha] = 0.25;
e = 1.602*10^-19;
kb = 1.38*10^-23;
T = 300;
d = 20*10^-6;
Gx = 10^17;
Dif = 4*10^-9;
Kr = 3.1*10^-9;
n0 = 10^21;
sol := {V, 
  D[#, x] &@
   NDSolve[{eq, bcs}, n[x], {x, 0, d}, MaxSteps -> 50000][[1, 1, 2]]}
res = Table[sol, {V, 0, 0.7, 0.05}];
ListPlot[Table[res, {x, 0, d}], Joined -> True, Frame -> True, 
 ImageSize -> 500]

Gx is a complicated function of x because I integrate a list as input. So, my result of DSolve is very big.

I need to stop the calculation of gradient as soon as it becomes zero. When the plot of n'[x] at x=0 cross the x axis (with values of V on it) my calculation should stop and the appropriate result (graph) should be displayed. I don't need the negative values of the derivative, so there is no need for Mathematica for further calculations. What is the best way to execute this?

But I have another problem there because my result is to high and it never changes the sign. In theory it should..

My result should look like:

Because I cannot post the images, you can see the desired graph on the link: https://www.google.si/search?hl=sl&rlz=1C1TEUA_enSI493SI493&q=current+voltage+characteristics+DSSC&bav=on.2,or.r_qf.&biw=1600&bih=785&um=1&ie=UTF-8&tbm=isch&source=og&sa=N&tab=wi&ei=YspCUb3GFYvUsgaTnYGoCA#imgrc=wX5CgAQmm0sdkM%3A%3BJbeOH9F0mOvjNM%3Bhttp%253A%252F%252Fars.els-cdn.com%252Fcontent%252Fimage%252F1-s2.0-S0025540812004047-gr9.jpg%3Bhttp%253A%252F%252Fwww.sciencedirect.com%252Fscience%252Farticle%252Fpii%252FS0025540812004047%3B384%3B314

On x axis is voltage, on y axis is current density, which I should calculate from the gradient at a point x=0: J=(dn[x]/dx)x=0. The current is related to a density gradient, the diffusion coefficient being the proportionality constant. But if I multiply the resulted derivative with Dif, it doesn't change enough. If I multiply the whole equation with e, which gives me correct units, the result is too small.

From the picture above is evident, that my values on y axis are changing to slowly.

Thank you for your answers!

Improve this question
$\endgroup$
11
  • $\begingroup$ Why do you say it is too high? $\endgroup$
    – wxffles
    Mar 13, 2013 at 21:45
  • 1
    $\begingroup$ Your code seems right. If your result for the derivative are much higher than you expected, perhaps you could check those constants. $\endgroup$
    – Dr. belisarius
    Mar 13, 2013 at 23:13
  • $\begingroup$ You can see analytically from the DSolve solution that this is correct. Why do you think that this is the wrong answer? $\endgroup$
    – Jonathan Shock
    Mar 14, 2013 at 6:29
  • $\begingroup$ I want to observe a current density. It is the gradient that I am calculating. But I should multiply my result by a proportionality factor Dif, which changes the result but not enough. The current density cannot have so high values. It is true, that also the units are not correct. I should multiply the result for gradient by elementary charge to get the current. But even I do this, the values are not ok. I must also find a way to stop the calculation, when my gradient (current density) is zero. I think If loop will do it. But first, I want to be sure, my calculation is correct.. $\endgroup$
    – Luka
    Mar 14, 2013 at 6:57
  • $\begingroup$ You can use the command AspectRatio if you just want to change the scale between the axes, for more see here $\endgroup$
    – Spawn1701D
    Mar 14, 2013 at 19:09

1 Answer 1

Reset to default
0
$\begingroup$

Please try the following 

eq = Gx + D[Dif*n'[x], x] -  Kr (n[x] - n0) == 0
bcs = {n'[d] == 0, n[0] == n0*Exp[(e α V)/(kb T)]};
α = 0.25;
e = 1.602*10^-19;
kb = 1.38*10^-23;
T = 300;
d = 20*10^-6;
Gx = 10^17;
Dif = 4*10^-9;
Kr = 3.1*10^-9;
n0 = 10^21;
sol := {V, 
D[#, x] &@NDSolve[{eq, bcs}, n[x], {x, 0, d}, MaxSteps -> 50000][[1, 1, 2]]}
res = Table[sol, {V, 0, 0.7, 0.05}];
Needs["PlotLegends`"]
ListPlot[(TakeWhile[#, #[[2]] > 0&]&/@(Table[#, {x, 0, d, d/100}]&/@Thread[{x, res[[All, 2]]}]))/.{}->Sequence[], Joined -> True, Frame -> True, ImageSize -> 500, PlotLegend -> res[[All, 1]], LegendLabel -> "V", AxesLabel -> {"x", "y"}]
Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.