Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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!

share|improve this question
    
Why do you say it is too high? –  wxffles Mar 13 '13 at 21:45
1  
Your code seems right. If your result for the derivative are much higher than you expected, perhaps you could check those constants. –  belisarius Mar 13 '13 at 23:13
    
You can see analytically from the DSolve solution that this is correct. Why do you think that this is the wrong answer? –  Jonathan Shock Mar 14 '13 at 6:29
    
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.. –  Luka Mar 14 '13 at 6:57
    
You can use the command AspectRatio if you just want to change the scale between the axes, for more see here –  Spawn1701D Mar 14 '13 at 19:09
show 6 more comments

1 Answer

up vote 0 down vote accepted

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 \[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}];
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"}]
share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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