Trying to solve an ODE (nonlinear, first-order)

[sjung-2]

Hello,

I just started using MATHCAD prime 3.0 without any experience with previous versions of MATHCAD. My first task is to solve an ODE (nonlinear, first-order: as analyzed by mathematica):

(dV(x)/dx)^2 - (dV(x)/dx|x=0)^2 = A*(exp(B*V(x))-1), with only two boundary conditions V(0) = 0 and V(d) = -Vd

For the details, please see the file attached below.

I guess the problem occurs because of the (dV(x)/dx|x=0)^2 term in the left-hand side of the ODE. But I do not have a clue how to fix this problem.

I tried to be as clear as possible.

Thank you very much,

Sungyeop

[Werner_E]

A first order ODE with 2 boundery conditions???? And Mathematica didn't complain?

You also have a unit problem: When calculating exp(X), X must be unitless, that means that B*V(x) is unitless and so the unit of V(x) will have to be Joule (1/UnitsOf(B)). The unit of the RHS is the the unit of A, which is Pascal. When V(x) are Joule, the LHS has the Unit Newton (J/m). You have different units on both sides which is not possible. And to top it your initial conditions say that V(x) should have the unit volt - the second actually says it has unit volt squared - you would have to delete the extra V you added as Vd already are volt. But you can't have two BC for an ODE of first order anyway.

The last argument of odesolve() should have the unit meter.

Some additional remarks:

Use the prime symbol instead of the differentiation operator.

Define the range for x (which is only meaningful to be used for the plot) after the solve solve block. To define a range with units, you have to give the unit for the first and second value as well.

According to Prime's help the ODE must be linear in the highest derivative for the numeric solver so do its job. This is not the case, so odesolve() will probably fail anyway.

[AndyWesterman]

From what I can reconstruct,I'm guessing that it is a semiconductor calculation.

first issue is with the "constants" that you have in the tables;

I assume that K[B is Boltxmann's constant...

If you go to the built-in constants in the ribbon at the top of the screen this is pre-defined as 'k'

From this I would question the values & units that you are using.

Feeding these changes to the equation for A then requires change to the definitions for E[b & σ, Both become voltages.

The units here still 'feel' wrong & I don't have enough information to make a guess, have a look at Nv

Because it is a conversion from P3 to P2 , I don't have any of the text fields , can you post a pdf of the sheet so that I can get a better understanding

Regards

Andy

[Werner_E]

Because it is a conversion from P3 to P2 , I don't have any of the text fields , can you post a pdf of the sheet so that I can get a better understanding

Here is the pdf.

How did you manage to make the P3 file at least partially readable in P2?

[AndyWesterman]

Hi Werner,

Thanks for the PDF, I think that my guesses are more or less right as far as I could go.

(PS: I'll outline the P3 -> P2 conversion in a seperate post - its messy!)

The definition for A is the next sticking point & that will be dependant on the derivation of the ODE later...

If V(x) is a voltage then

d/dx( V(x) ) is in V/m & the ODE would require A to have dimensions of V^2 / m^2 ...

Sungyeop, can you advise what the ODE is modeling & where the equation was derived?

Regards

Andy

[sjung-2]

Here, I converted units into [J] rather than using [eV].

Would it be helpful, Westerman?

The equation is Poisson's equation: d^2V(x)/dx^2 = - q*p(x)/eps, ....... ODE(1)

where p(x) = -(q*Nv)/eps * exp[-1/(kT)*{Eb - sigma^2/(2*k*T)}] * exp(-q*V(x)/(k*T))

V(x) is electrical potential [V]

x: physical dimension [m]

q: elemental charge [C]

Nv: density of states [cm^(-3)]

Eb: energy barrier [eV or J]

sigma: disorder factor [eV or J]

eps: permittivity [F*m^(-1)]

By setting F(x) = - dV(x)/dx, ODE(1) becomes,

dV(x)/dx = [A*exp(B*V(x))+V'(0)^2]^(1/2) ........ ODE(2)

where

A = - 2 * (kT)/eps * Nv * exp[-1/(kT)*{Eb-sigma^2/(2*k*T)}]

B = - q/(k*T)

Thus, A would have [J*K-1]*[K]*[F-1*m]*[m-3] = [J*F-1*m-2] = [V2m-2]

Having said that, I need to solve either ODE(1) or ODE(2). The problem is, for both cases, the fact that I do not know the value of V'(0). Knowing only V(0) = 0, V(d) = -0.6.

Regards,

Sungyeop

[Werner_E · Accepted solution]

Now the units are consistent. You have some minor errors as you tend to add the unit to a variable which already is assigned a unit (e.g. use just "d", not "d m"). Also for the prime symbol you are not allowed to just type an apostrophe, you have to use the symbol from the menu or use the keyboard shortcut.

I added the minerr solve block Alan was talking about to provide different values of C (=V'(d)) and try to chose C so that V(d)=-V.d.

Its no surprise that we cannot get negative values for V(d) (you demand -0.6V), no matter how we would chose C. The ODE is setup that V'(x) equals a square root which is positive and so the function is ascending. Chosing the ODE to V'(x)= minus square root... makes the odesolveblock fail for reasons unknown to me.

There is also a further problem as the values for the endpoint V(d) seems not to change continually when C is changed but apruptly in steps. That make the minerr() solveblock to just return the guess value, no matter if you demand Vd(C)=-V.d or Vd(C)=+V.d. I am not sure if thats because of a to small numer of steps for the odesolve - I guess, thats not the reason. Increasing the number of steps for odesolve does not change that behaviour up to the point where the minerr() block would fail.

I add my attempt anyway - maybe someone else want to investigate some more time.

[AndyWesterman]

The part that I'm struggling with at the moment is the translation from

to

I've followed the equations as stated from the last e-mail & at the mon=ment I get numbers >10^307 so obviously something wrong.

regards

Andy

[Werner_E]

I had not looked at the mail in detail, but I gues Sungyeop would have to explain his substitution in more detail. I am not sure if his intitial second order ODe1 is correct.

[AlanStevens]

When you've fixed the problems Werner has noted you should move the dV/dx|x=0 term (which is a constant) to the right hand side and take the square root of both sides. You will then have a first-order ODE. Use V(0) as the initial condition. Presumably, you don't know the value of dV/dx|x=0, which is why the problem is formulated as it is. In this case you will have to guess values for dV/dx|x=0 and solve repeatedly until you find one that results in V(d) = -Vd.

Set up the ODE as a function with dVdx|x=0 as an input and call it from another solve block that minimizes the error between V(d) and -Vd.

(Can't help with an explicit worksheet as I don't have Prime 3)

Alan

[sjung-2]

Thank you very much,

It becomes much more clear.

Can you please let me how to set up an iterative resolution in MATHCAD? With some reference dealing with it?

Regards,

Sungyeop

[AlanStevens]

Go to the help files and search for minerr. There should be an example of how to parameterize a solve block there.

Alan

[sjung-2]

After a long absent, I would like to make some further notes. Thank you very much in advace.

1) Negative boundary condition V(d) = -0.6 V

As Werner mentioned, mathematically, V(d) cannot take a negative value. And, I have noticed that a negative V(d) is also physically non-sense. Sorry for my carelessness. For it to be physically and mathematically meaningful, V(d) should always be 0.6 V. In fact, in the problem, V(0) is the electrical potential imposed by Au(gold) electrode and V(d) is that of Al(Aluminium).

2) The slow change of V(d) in response to the change of C(the electric field, which is dV/dx).

As mentioned in the last part of 1), we are dealing with a problem where the electrical potential at two electrodes V(0) and V(d) are fixed at 0 and 0.6V, respectively. Thus, in principle, V(0) and V(C) should not be changing in response to the change of C(the electric field).

The ODE solver was written with one boundary condition. It would be better if we can impose two boundary conditions V(d) = 0.6 V as well as V(0) = 0 V. But I have no clue, for now, how to do that.

Anyway, the proposed solution was very useful.

The variable that I am trying to sweep is 'sigma' from 0.05 to a value where V(d) restes on 0.6 V, while struggling to extract data from the plot using excel component. (recently uploaded an other thread on that problem. If time allows, please have a look.)

I hope this response would help you understand the problem better. Thank you.