# Self-Consistent solution of Boltzmann eqn. using Matlab

I am trying to solve Boltzmann eqn:

$\frac{e}{m}\tau\vec{E}_x\cdot \nabla_\vec{v}{f(\vec{v})}=f(\vec{v})-f_0(\vec{v})$

where

$f_0(v)=\exp{\frac{-m \vec{v}^2}{T}}$ and $\vec{v}=\{\vec{v}_x,\vec{v}_y\}$ (i.e., $\vec{v}$ is in 2D.)

For simplicity, let's assume $\vec{v}_y=0$. Assume $\tau$ and $T$ are constants. Then for any given $\vec{v}$ and $|E|$, $f_0(\vec{v})$ is completely known.

I am trying to set up a self-consistent code in Matlab to solve this equation for any given set of $\vec{v}$. However, I am not sure how to get Matlab to solve this? the ODE command, ode45, is not exactly for this task. Any ideas?

EDIT: I need to do a simulation, so I will choose the values of $\vec{v}$ and do the simulation accordingly. All other constants are assumed known ($e,m,\tau,T$ and $|E|$).

• Welcome to SciComp! As posed, your question is unclear. Is everything given besides $v$? It looks like you'll need to solve a first order steady-state partial differential equation, in which case, you'll have to discretize your equation in one direction before using ode45. It's not clear that doing so is the best course of action; you might be better off discretizing both directions instead and solving the resulting nonlinear equations. I can't say what type of discretization would work best without knowing more about the problem. – Geoff Oxberry May 6 '14 at 1:54
• Thanks, I have edited my post. Yes everything is known, and the values of $\vec{v}$ is chosen by us to do a simulation. What is remaining is to do the self-consistent solution. – student1 May 6 '14 at 1:58

To give background of the problem the general Boltzmann Transport Equation for electrons in an electric field is

$\frac{\partial f}{\partial t} + \vec{v} \cdot \nabla_{\vec{r}} f - \frac{e}{m} \vec{E} \cdot \nabla_{\vec{v}}f = \left(\frac{\partial f}{\partial t}\right)_{\mathrm{coll}}$

where $f$ is the phase space velocity distribution function. $f_0(\vec{v})$ is the steady state solution with no field (i.e. a maxwellian velocity distribution).

Then a (usually small) field is applied in one direction ($E_x$) which perturbs the velocity distribution of the electrons. At this point, the distribution changes slightly, and if we apply steady state $\left(\frac{\partial f}{\partial t} = 0\right)$and isotropic $\left(\nabla_{\vec{r}} f = 0\right)$ we are left with:

$\frac{e}{m} \vec{E} \cdot \nabla_{\vec{v}}f = \left(\frac{\partial f}{\partial t}\right)_{\mathrm{coll}}$

now we simply assume that $\left(\frac{\partial f}{\partial t}\right)_{\mathrm{coll}}$ is proportional to the perturbation from equilibruim, $f(\vec{v}) - f_0(\vec{v})$, with a proportionality constant of $\frac{1}{\tau}$. Defining the $x$ axis along the field gives the equation from the question:

$\frac{e}{m} \tau |E| \frac{d f(\vec{v})}{d v_x} = f(\vec{v}) - f_0(\vec{v})$

Now due to the steady state being isotropic and us defining the x direction along the field we get:

$\frac{e}{m} \tau |E| \frac{d f(v_x)}{d v_x} = f(v_x) - f_0(v_x)$

This form of the equation has an exact analytic solution (it is just a 1st order ODE). Let

$a = \frac{e}{m} \tau |E|$

and the ODE becomes

$f'(v_x) - \frac{1}{a} f(v_x) = - \frac{1}{a} f_0(v_x)$

which has a general solution:

$f(v_x) = C \exp(\frac{v_x}{a}) +\exp(\frac{v_x}{a}) \int_1^x - \frac{1}{a} \exp{\frac{\xi}{a}} f_0(\xi) \mathrm{d}\xi$

(shameless pulled from wolfram alpha at the moment)

I don't have my notes with me at the moment to make sure I have done this all exactly correctly, but I will update this tomorrow when I get to the office. Im also fairly certain that the final equation simplifies but I can't remember how.

EDIT:

So I went back through my notes, and the problem I remembered from stat thermo had a few other small assumptions that don't apply in your question (it converted to energy as well). Using wolframalpha and unity for the constants, I got an exact solution, although I am not happy with it (the inverse error function doesn't seem appropriate).

Even still, it should be possible to solve the above equation numerically as you originally asked using a standard finite difference formulation. The boundary conditions should be that $f$ goes to 0 at both positive and negative infinity and $\int_{-\infty}^\infty f dv_x = 1$.

• I happen to use a slightly more in depth set of assumptions for this exact equation in my research, so it was more for me practicing than anything. – Godric Seer May 6 '14 at 3:11
• The solution to the 1st order ODE will involve the error function, I wonder whether this is reasonable. – student1 May 6 '14 at 3:47
• I added the analytic solution. In about 12 hours I can go through my old notes and sort out where exactly this analysis ends up, but I am 95% certain that the point of assuming the collision term is proportional to the perturbation is that it can be solved exactly. – Godric Seer May 6 '14 at 3:53
• @student1 I have edited my post. Apparently I slightly misremembered my previous problem and the assumptions it made. – Godric Seer May 6 '14 at 15:26
• Thanks. Yes, erfc(x) does not look appropriate here to me either. The problem in question clearly asks for a solution using self-consistency, for various electric fields (i.e. E is a variable, and so is v, so we get 3D plot). – student1 May 6 '14 at 21:04