User:Egm6322.s09.Three.nav

PROJECT # 2
Project

The word linear comes from the latin word ‘linearis’ which means created by lines.

'''Proof to show that the order of double differentiation is independent of the variables i.e uxy= uyx Assume that there exist uniform number of intervals of length 'h' in both 'x' and 'y' directions'''

Then

ux= $$\lim_{h\to 0}\frac {u(x+h,y)-u(x,y)}{h}$$

uy= $$\lim_{h\to 0}\frac {u(x,y+h)-u(x,y)}{h}$$

uxy= $$ \lim_{h\to 0}\frac {\frac{\partial u(x+h,y) }{\partial y}-\frac{\partial u(x,y)}{\partial y}}{h}$$

$$\Rightarrow \lim_{h\to 0}\frac{\left [\frac{u(x+h,y+h)-u(x+h,y)}{h} \right ]-\left [\frac{u(x,y+h)-u(x,y)}{h} \right ]}{h}$$

$$\Rightarrow \lim_{h\to 0}\frac{u(x+h,y+h)-u(x+h,y)-u(x,y+h)+u(x,y)}{h^{2}}$$---(1) uyx= $$ \lim_{h\to 0}\frac {\frac{\partial u(x,y+h) }{\partial x}-\frac{\partial u(x,y)}{\partial x}}{h}$$

$$\Rightarrow \lim_{h\to 0}\frac{\left [\frac{u(x+h,y+h)-u(x,y+h)}{h} \right ]-\left [\frac{u(x+h,y)-u(x,y)}{h} \right ]}{h}$$

$$\Rightarrow \lim_{h\to 0}\frac{u(x+h,y+h)-u(x+h,y)-u(x,y+h)+u(x,y)}{h^{2}}$$---(2) Hence As is evident from (1) and (2), the order of double differentiation is indeed independent of the variables.

PROJECT#3
To do: Transform the Laplace Equation into Polar coordinates and compare with published results

The Laplace's Equation is often encountered in electrostatics, heat and mass transfer theory and in other areas of mechanics and physics. Simply worded, the Laplace's Equation states that div(grad u)=0.

It has the following mathematical form:

In Cartesian coordinates (x,y)$$: \ \ u_{xx}+u_{yy}= 0 $$

In Polar coordinates (r,$$\theta$$)$$: \ \ \ \ \ u_{rr}+\frac{u_{r}}{r}+\frac{u_{\theta\theta} }{r^{2}}=0 $$

where x= rcos($$\theta$$) and y= rsin($$\theta$$)

$$ Recognize\ that\ u_{xx}+u_{yy}= \begin{bmatrix} \partial _{x}\ \partial_{y} \end{bmatrix} \begin{bmatrix} \partial _{x} \\ \partial_{y}\end{bmatrix} \left \{u \right \} =\begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}JJ^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix}\left \{u \right \}$$

where $$ J= \begin{bmatrix} \frac {\partial r}{\partial x} \ \frac {\partial r}{\partial y}\\ \frac {\partial \theta }{\partial x} \ \frac{\partial \theta}{\partial y} \end{bmatrix}$$

Denote cos($$\theta$$)=C and sin($$\theta$$)=S

$$ \Rightarrow J= \begin{bmatrix} \ C \ S \\\frac{-S}{r} \ \frac{C}{r} \end{bmatrix}$$

Differentiating by parts

$$\begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}JJ^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix}\left \{u \right \}= \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}\left \langle JJ^{T}\right \rangle\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}+\left \langle \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix} J\right \rangle J^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}$$---(A)

where the terms in the angled brackets are kept fixed throughout differentiation.

Evaluating the RHS of equation (A) term by term '''First term= $$\begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}\left \langle JJ^{T}\right \rangle\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \} $$(1)

$$ JJ_{T}=\begin{bmatrix} \ C \ S \\ \frac{-S}{r} \ \frac{C}{r} \end{bmatrix}\begin{bmatrix} \ C \ -S/r \\ S \ \ C/r \end{bmatrix}= \begin{bmatrix} 1 \  0 \\  0 \ \frac{1}{r^{2}} \end{bmatrix}$$

$$\Rightarrow \ (1)= \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}\left \langle \begin{bmatrix} 1 \ \ \ 0 \\ \ 0 \ \frac{1}{r^{2}} \end{bmatrix} \right \rangle\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}$$

$$\Rightarrow\ (1)=\begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}\begin{bmatrix} \partial_{r} \\ \frac{\partial_{\theta}}{r^2}\end{bmatrix}\left \{u \right \}$$

$$\Rightarrow\ (1)= \frac {\partial^2{u}}{\partial{r^2}}+\frac{1}{r^2}\frac{\partial^2{u}}{\partial{\theta^2}} $$ ''' Second term= $$\left \langle \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix} J\right \rangle J^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}'''$$-(2)

By definition of the Jacobian, note that the term in the angled brackets is the same as $$\begin{bmatrix} \partial_{x} \ \partial_{y} \end{bmatrix}$$

$$\Rightarrow\ (2)= \left \langle \begin{bmatrix} \partial_{x} \ \partial_{y} \end{bmatrix}J^{T} \right \rangle \begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left\{u\right \}$$

Evaluating the term in angled brackets for equation(2)

$$\left \langle \begin{bmatrix} \partial_{x} \ \partial_{y} \end{bmatrix}J^{T} \right \rangle= \begin{bmatrix} \partial_{x}C+\partial_{y}S \ -\partial_{x}\frac {S}{r}+\partial_{y}\frac{C}{r} \end{bmatrix}$$

$$= \begin{bmatrix} \frac{1}{r} \ 0 \end{bmatrix}$$

$$\Rightarrow\ (2)= \begin{bmatrix} \frac{1}{r} \ 0 \end{bmatrix}\begin{bmatrix} \partial_{r} \\ \partial_{\theta} \end{bmatrix} \left \{u \right \} $$

$$\Rightarrow\ (2)= \frac{1}{r} \partial_{r}u$$

Finally putting terms (1) and (2) together into equation (A)

$$u_{xx}+u_{yy}= \begin{bmatrix} \partial _{x}\ \partial_{y} \end{bmatrix} \begin{bmatrix} \partial _{x} \\ \partial_{y}\end{bmatrix} \left \{u \right \} =\begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}JJ^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix}\left \{u \right \}$$

$$\Rightarrow\ (A)= \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}JJ^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix}\left \{u \right \}= \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix}\left \langle JJ^{T}\right \rangle\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}+\left \langle \begin{bmatrix} \partial_{r}\ \partial_{\theta} \end{bmatrix} J\right \rangle J^{T}\begin{bmatrix} \partial_{r}\\ \partial_{\theta} \end{bmatrix} \left \{u \right \}$$

$$\Rightarrow\ u_{xx}+u_{yy}= \frac {\partial^2{u}}{\partial{r^2}}+\frac{1}{r^2}\frac{\partial^2{u}}{\partial{\theta^2}}+\frac{1}{r}\frac {\partial{u}}{\partial{r}}$$

Hence we have proved, using transformation of coordinates, that:

$$ u_{xx}+u_{yy} =\ u_{rr}+\frac{1}{r}u_{r}+\frac{1}{r^2}u_{\theta\theta}$$

Comparing with published results, we see that our derived result matches perfectly.

Rotation of Coordinates
What is the form of the Jacobian when implementing a rotation of coordinates?

Consider one coordinate system, x-y (cartesian) and another coordinate system, $$\bar{x}-\bar{y}$$, obtained by rotating x-y clockwise by angle $$\theta$$ (as shown in the figure)

Assuming $$\widehat{\underline{i}}$$ and $$\widehat{\underline{j}}$$ are the unit vectors in the new coordinate system, $$\hat{i}$$ and $$\hat{j}$$ are the unit vectors in the old coordinate system

Then

$$ \widehat{\underline{i}}= icos\left(\theta \right)+ jsin\left(\theta \right)$$

$$ \widehat{\underline{j}}= -isin\left(\theta \right)+ jcos\left(\theta \right)$$

Comparing the above vector equations with those given in the example of Affine Mapping  considered previously,

$$\begin{matrix} \widehat{\underline{i}}=m \widehat{i}+n \widehat{j}\\ \widehat{\underline{j}}=p \widehat{i}+q \widehat{j} \end{matrix}$$

We see that m= cos$$\theta$$, n= sin$$\theta$$, p= -sin$$\theta$$, q= cos$$\theta$$

Hence the Jacobian ( J ) is given by

$$\underline{J}= \underline{E}= \begin{bmatrix} \frac{q}{\left (mq-np \right )}&\frac{-p}{\left (mq-np \right )} \\ \frac{m}{\left (mq-np \right )}&\frac{-n}{\left (mq-np \right )} \end{bmatrix}$$

$$\because$$mq-np= $$ \left(cos\theta \times cos\theta \right)$$ - $$\left(sin\theta \times -sin\theta \right)$$ = $$ cos^{2}\theta+ sin^{2}\theta$$ = 1

$$ \therefore \underline{J}= \underline{E}= \frac {1} {mq-np}\begin{bmatrix} \ \ q \ \ -p \\ -n \ \ \ m \end{bmatrix} = \frac {1} {1}\begin{bmatrix} cos\theta \ \ \ sin\theta \\ -sin\theta \ \ cos\theta \end{bmatrix} = \begin{bmatrix} cos\theta \ \ \ sin\theta \\ -sin\theta \ \ cos\theta \end{bmatrix}$$

PROJECT#5
Particularizing the general solution:

Case (i): Domain $$ \mho $$ is a full circle i.e. $$\theta \ \epsilon \ [0,2\pi]$$

$$ \Rightarrow$$ Solution has to be periodic w.r.t $$\theta$$

$$ \Rightarrow u(r,\theta+2k\pi, t) = u(r,\theta,t) $$, for any non-negative integer 'k'.

This implies that solution for $$\theta$$($$\theta$$), which consists of sines and cosines (periodic functions) is only possible only for integers n= 1,2,...

$$ \Rightarrow \sqrt{\rho}= n $$

$$ \Rightarrow \rho= n^{2} $$

Also replace the continuous integral $$ \int d\rho$$ by a discrete sum $$\sum_{0}^{n}$$ in Eqn.(2)

Now the equation becomes

$$ u(r,\theta,t)= \int_{0}^{\infty} d\lambda \left [\sum_{0}^{n} exp(-\lambda t)\left (\bar{B}(\rho, \lambda)sin(\sqrt{\rho}\theta) + \bar{C}(\rho,\lambda)cos(\sqrt{\rho}\theta) \right)\left (D(\rho, \lambda)J_{\sqrt{\rho}}(\sqrt{\lambda}r)+(E(\rho, \lambda)Y_{\sqrt{\rho}}(\sqrt{\lambda}r \right) \right] $$

An important point to remember here is that the above solution, obtained as a result of assuming periodicity, is applicable to the annulus domain.

Q: What if the domain is a solid disk, instead of an annulus?

In the plots of $$ J_{n}(z)$$ and $$Y_{n}(z)$$ Vs z, it was seen that $$Y_{n}(z) \to -\infty \ as\ z\ \to \ z^{+} $$

For the solution to remain meaningful, this cannot be allowed to happen. Hence we force $$E(\rho, \lambda)$$ = 0.

Hence the equation reduces to

$$ u(r,\theta,t)= \int_{0}^{\infty} d\lambda \left [\sum_{0}^{n} exp(-\lambda t)\left (\bar{\bar{B}}(\rho, \lambda)sin(\sqrt{\rho}\theta) + \bar{\bar{C}}(\rho,\lambda)cos(\sqrt{\rho}\theta) \right)J_{\sqrt{\rho}}(\sqrt{\lambda}r)\right] $$

where $$\bar{\bar{B}}= \bar{B}(\rho,\lambda)D(\rho,\lambda)$$ and

$$\bar{\bar{C}}= \bar{C}(\rho,\lambda)D(\rho,\lambda)$$

Remembering that $$\sqrt{\rho} = n$$, the above equation can be rewritten as

$$ u(r,\theta,t)= \int_{0}^{\infty} d\lambda \left [\sum_{0}^{n} exp(-\lambda t)\left (\bar{\bar{B}}(\rho, \lambda)sin(n\theta) + \bar{\bar{C}}(\rho,\lambda)cos(n\theta) \right)J_{n}(\sqrt{\lambda}r)\right] $$

DERIVATION OF FOURIER COEFFICIENTS
How do we derive the fourier coeffcients C0, An and Bn?

An orthogonal basis: {1, cosm$$\theta$$, sinm$$\theta$$} is used.

As defined above, two functions f($$\theta$$), g($$\theta$$) are orthogonal if $$\int_{0}^{2\pi}f(\theta).g(\theta) d\theta= 0$$

$$\blacktriangleright$$Eg. Consider f($$\theta$$)= cos ($$m\theta$$) and g($$\theta$$)= cos($$n\theta$$)

$$ \int_{0}^{2\pi} cos(m\theta).cos(n\theta)d\theta= \int_{0}^{2\pi} \frac{1}{2}\left (cos \left ((m+n)\theta \right)+cos \left((m-n)\theta \right) \right)d\theta $$

$$ = \left [ \left (\frac{1}{2(m+n)} (sin \left ((m+n)\theta \right) \right)+ \left (\frac{1}{2(m-n)} sin \left((m-n)\theta \right) \right) \right]_{0}^{2\pi}$$

$$ = \begin{Bmatrix} 0, if\ m\ \neq n\\ 2\pi, if\ m\ = n

\end{Bmatrix}$$

For more details, here is a link that explains the math in greater detail. Similarly using '1' as the basis function and simplifying, we get $$\int_{0}^{2\pi}cos\theta\left \{1 \right \} d\theta \ or\ \int_{0}^{2\pi}sin\theta \left \{1 \right \} d\theta = 0 $$

Using this knowledge, we proceed to determine the fourier coefficients C0, An and Bn, using the orthogonal basis {1, cosm$$\theta$$, sinm$$\theta$$}

Consider the equation $$T(r,\theta)=C_{0}+\sum_{n=1}^{\infty }r^{n} \left \{A_{n}cosn\Theta +B_{n}sinn\Theta \right\}$$-(1)

Say we have the Boundary condition T(r=a, $$\theta$$)= ($$ T^{*}\left(\theta \right)$$)

Then, $$T^{*}\left(\theta \right) = C_{0}+\sum_{n=1}^{\infty }a^{n} \left \{A_{n}cosn\Theta +B_{n}sinn\Theta \right\}$$(2)

$$\blacktriangleright$$At n=0, $$T^{*}\left(\theta \right) = C_{0}$$

Multiplying through by {1} and integrating over [0,2$$\pi$$]

$$\Rightarrow \int_{0}^{2\pi}T^{*}\left(\theta \right) d\theta = \int_{0}^{2\pi} C_{0} d\theta$$

$$\Rightarrow \int_{0}^{2\pi}T^{*}\left(\theta \right) d\theta = C_{0} \left (\theta \right)_{0}^{2\pi}$$

$$\Rightarrow C_{0}= \frac {1}{2\pi} \int_{0}^{2\pi}T^{*}\left(\theta \right) d\theta $$

$$\blacktriangleright$$Solving for An in Equation (2), multiply through by cos(m$$\theta$$) and integrate over [0,2$$\pi$$].

Evaluating each term in the resultant equation,

Term on the LHS= \int_{0}^{2\pi}T^{*}\left(\theta \right)cosm\theta d\theta

First term on RHS= $$\int_{0}^{2\pi}C_{0}cosm\theta d\theta $$

$$= C_{0} \int_{0}^{2\pi} \left \{1 \right \} cos m\theta d\theta $$

= 0, by definition of orthogonality

Second term on RHS= $$ \int_{0}^{2\pi} \sum_{n=1}^{\infty}a^{n}A_{n}cosn\theta cosm\theta d\theta $$

$$ = \sum_{n=1}^{\infty} \int_{0}^{2\pi}a^{n}A_{n}cosn\theta cosm\theta d\theta $$

$$ = \begin{Bmatrix} 0, if\ m\ \neq n\\ a^{n}\times A_{n} \times 2\pi, if\ m\ = n

\end{Bmatrix}$$ (by definition of orthogonality)

$$\Rightarrow \int_{0}^{2\pi} \sum_{n=1}^{\infty}a^{n}A_{n}cosn\theta cosm\theta d\theta = a^{n}\times A_{n} \times 2\pi $$

Hence it is seen that though it was assumed m $$\epsilon$$ $$ \Re$$, simplification of the second term in (2) determines that m=n.

Third term on RHS= $$ \int_{0}^{2\pi} \sum_{n=1}^{\infty}a^{n}A_{n}sinn\theta cosm\theta d\theta $$

= 0, by definition of orthogonality

$$\therefore $$ (2) multiplied through by cos($$m\theta$$) or cos($$n\theta$$) and integrated over [0,2$$\pi$$] reduces to

$$\int_{0}^{2\pi}T^{*}\left(\theta \right)cosn\theta d\theta = a^{n}\times A_{n} \times 2\pi$$

$$\Rightarrow A_{n}= \frac {1}{a^{n}\times 2\pi} \int_{0}^{2\pi}T^{*}\left(\theta \right)cosn\theta d\theta $$

We solve similarly for Bn in (2). Multipling through with sin($$m\theta$$) and integrating over [0,2$$\pi$$], it is seen that the first and second terms = 0, by definition. Third term reduces to $$b^{n}\times B_{n} \times 2\pi$$. Then the modified Equation(2) becomes

$$\int_{0}^{2\pi}T^{*}\left(\theta \right)sin n\theta d\theta = a^{n}\times B_{n} \times 2\pi$$

$$\Rightarrow B_{n}= \frac {1}{b^{n}\times 2\pi} \int_{0}^{2\pi}T^{*}\left(\theta \right)sin n\theta d\theta $$

PROJECT#6
The Leibniz Integral Rule, named after Gottfried Leibniz accomplishes differentiation of an integral. It tells us that

$$\frac{\mathrm{d} }{\mathrm{d} x}\left(\int^{B(x)}_{A(x)}F(x,\xi)d\xi \right)= \int^{B(x)}_{A(x)}\frac{\partial F(x,\xi)}{\partial x}d\xi+\left(F(x, \xi= B(x))\frac{\mathrm{d} B}{\mathrm{d} x} \right)- \left(F(x, \xi= A(x))\frac{\mathrm{d} A}{\mathrm{d} x} \right)$$ --(1)

This can be derived from the fundamental theorem of calculus as shown in this Wikipedia link. Leibniz_integral_rule

The Leibniz Theorem is applied as a tool in deriving most of the conservation laws in Fluid Mechanics. Consider a control volume that is neither fixed nor moving, but is simply bounded by surfaces a(t) and b(t). These surfaces are moving with a velocity different from that of the local fluid velocity. Then the differential of the integral, as given by the Leibniz theorem, is written as

$$\frac{\mathrm{d} }{\mathrm{d} t}\left(\int^{b(t)}_{a(t)}F(x,t)dx \right)= \int^{B(x)}_{A(x)}\frac{\partial F(x,t)}{\partial t}dx+\left(F(x= b(t),t)\frac{\mathrm{d} B}{\mathrm{d} t} \right)- \left(F(x= a(t),t))\frac{\mathrm{d} A}{\mathrm{d} t} \right)$$-(2)

where a(t), b(t) form the boundaries of the control volume considered and are moving at a speed dA/dt and dB/dt respectively.

Multiplying through by dS (the surface area of a differential control volume) and generalizing further, we obtain the following form:

$$\frac{\mathrm{d} }{\mathrm{d} t}\left(\int^{b(t)}_{a(t)}F(x,t)dV \right)= \int^{b(t)}_{a(t)}\frac{\partial F(x,t)}{\partial t}dV+\int_{S(t)}\mathbf{dS.u_{s}}F$$--(3)

The surface integral over area S(t) in (3) includes the contributions from both the boundaries, thus excluding the need for the previous form of the equation (2).

For a control volume that is fixed, i.e. u= 0, (3) becomes

$$\frac{\mathrm{d} }{\mathrm{d} t}\left(\int^{b(t)}_{a(t)}F(x,t)dV \right)= \int^{b(t)}_{a(t)}\frac{\partial F(x,t)}{\partial t}dV$$--(4)

For a material volume i.e. a control volume in which the boundaries are moving at a velocity equal to the fluid velocity, u, (3) becomes

$$\frac{\mathrm{D} }{\mathrm{D} t}\left(\int^{b(t)}_{a(t)}F(x,t)dV \right)= \int^{b(t)}_{a(t)}\frac{\partial F(x,t)}{\partial t}dV+\int_{S(t)}\mathbf{dS.u_{s}}F$$--(5)

Eqn.(5) is sometimes referred to as the Reynold's Transport Theorem'. We use the term D/Dt to convey that we are defining the property of a material control volume. This theorem is often used in formulating the basic conservation laws in fluid mechanics.

Applying Gauss' theorem, which states that $$\oint_{S}\mathbf{X.dS}= \int_{V}(\nabla.X)dV$$, (5) can be rewritten as

$$\frac{\mathrm{D} }{\mathrm{D} t}\left(\int^{b(t)}_{a(t)}F(x,t)dV \right)= \int^{b(t)}_{a(t)}\frac{\partial F(x,t)}{\partial t}dV+\int_{V(t)}\nabla.(u_{s}F)dV$$(6)

Consider F as a physical property

F= $$\rho$$ gives us the equation for Conservation of Mass. For a steady flow system, the mass of the system i.e.$$\rho$$dV remains constant with time. Hence from(6)

$$0= \int^{b(t)}_{a(t)}\frac{\partial \rho}{\partial t}dV+\int_{V(t)}\nabla.(\rho u_{s})dV$$

(or)

$$\int_{V(t)}\left (\frac{\partial \rho}{\partial t}+\nabla.(\rho u) \right) dV = 0 $$---(7)

Similarly F= linear momentum= $$\rho$$v gives us the equation for Conservation of Momentum.

$$\int_{V(t)}\left (\frac{\partial \rho v}{\partial t}+\nabla.(\rho vu) \right) dV = 0 $$

Hence we see the relevance of the Leibniz Integral rule in the context of Fluid Mechanics and the fundamental role it plays in deriving the basic laws that the complex science of Fluid Mechanics is built on.

PROJECT#7
Problem:  Given the general,modified form of a PDE: $$\lambda_{1}\bar{x}^2+\lambda_{2}\bar{y}^2+d\bar{x}+e\bar{y}+f=0$$--(1)

Derive the canonical form representing a parabola from (1).

The canonical form representing a parabola is given by the following equation

$$\bar{\xi}^{2}- \bar{\eta} = 0$$ --(2)

where $$\bar{\xi}$$ and $$\bar{\eta}$$ represent coordinates in a different coordinate system.

So how do we derive (2) from (1)? We apply rotation once and translation of coordinates twice to $$\bar{x}$$ and $$\bar{y}$$.

Step(1)- First Translation

Define $$\begin{bmatrix}\bar{\bar{x}}\\ \bar{\bar{y}}\end{bmatrix}= \begin{bmatrix}\bar{x}+r \\ \bar{y}\end{bmatrix}$$

Substituting $$\bar{\bar{x}}$$ and $$\bar{\bar{y}}$$ in (1), we get

$$\lambda_{1}(\bar{\bar{x}}-r)^2+\lambda_{2}\bar{\bar{y}}^2+d(\bar{\bar{x}}-r)+e\bar{\bar{y}}+f=0$$

$$\Rightarrow \lambda_{1}\bar{\bar{x}}^{2}+\lambda_{2}\bar{\bar{y}}^2+\bar{\bar{x}}(d-2\lambda_{1}r)+e\bar{\bar{y}}+\lambda_{1}r^{2}+f=0$$--(3)

For equation (3) to resemble (2), we need the terms: $$\bar{\bar{y}}^2, \bar{\bar{x}} \ and \ \lambda_{1}r^{2}+f$$ to vanish.

Note that for a parabola, the equation of conics: $$Ax^2+Bxy+Cy^2+\bar{d}x+\bar{y}+\bar{f}=0$$, $$B^2-4AC=0$$ i.e only one eigenvalue exists. Thus assume $$\lambda_{2}$$= 0.

Also in order for the $$\bar{\bar{x}}$$ term to vanish, we need its coefficient to be equal to 0 $$\Rightarrow d- 2\lambda_{1}r=0$$

$$\Rightarrow r= \frac{d}{2\lambda_{1}r}$$

Substituting $$\lambda_{2}$$ and r in (3), the equation reduces to $$\lambda_{1}\bar{\bar{x}}^{2}+e\bar{\bar{y}}+\frac{d}{4\lambda_{1}^{2}}+f=0$$-(4)

Denote g= $$\frac{d}{4\lambda_{1}^{2}}+f$$. Thus (4) reduces to

$$\Rightarrow \lambda_{1}\bar{\bar{x}}^{2}+e\bar{\bar{y}}+g=0$$-(5)

Step(2)- Rotation

Eqn.(4) can be rewritten as $$\left \lfloor \bar{\bar{x}} \ \bar{\bar{y}} \right \rfloor\begin{Bmatrix} \lambda_{1} \ 0 \\ 0 \ \ 0 \end{Bmatrix}\begin{bmatrix}\bar{\bar{x}}\\ \bar{\bar{y}}\end{bmatrix}+\left \lfloor 0 \ e \right \rfloor\begin{bmatrix}\bar{\bar{x}}\\ \bar{\bar{y}}\end{bmatrix}+g=0$$ (6)

Now apply rotation of coordinates by defining a new coordinate system $$(\xi, \eta)$$ where $$\begin{bmatrix} \xi \\ \eta \end{bmatrix}= \mathbf{J_{\beta}}\begin{bmatrix}\bar{\bar{x}}\\ \bar{\bar{y}}\end{bmatrix}$$

where $$\mathbf{J_{\beta}}= \begin{bmatrix}\frac{1}{\sqrt{\lambda_{1}}} \ \ \ 0 \\ \ \ \ 0 \ \ \ \frac{1}{e}\end{bmatrix}$$

Due to $$\mathbf{J_{\beta}}$$ being an orthogonal matrix, $$\mathbf{J_{\beta}}^{-1}= \mathbf{J_{\beta}}^{T}= \mathbf{J_{\beta}}$$

$$\Rightarrow \begin{bmatrix}\bar{\bar{x}}\\\bar{\bar{y}} \end{bmatrix}= \mathbf{J_{\beta}}^{-1}\begin{bmatrix}\xi\\\eta \end{bmatrix}= \mathbf{J_{\beta}}\begin{bmatrix}\xi\\\eta \end{bmatrix}$$

Substituting this rotation of coordinates in (6), we get

$$\begin{bmatrix}\xi\ \eta \end{bmatrix} \mathbf{J_{\beta}} \begin{Bmatrix}\lambda_{1} \ 0 \\ 0 \ \ 0 \end{Bmatrix} \mathbf{J_{\beta}}\begin{bmatrix}\xi\\\eta \end{bmatrix}+\left \lfloor 0 \ e \right \rfloor \mathbf{J_{\beta}}\begin{bmatrix}\xi\\\eta \end{bmatrix}+g=0$$

$$\Rightarrow \begin{bmatrix}\xi\ \eta \end{bmatrix} \begin{bmatrix}\frac{1}{\sqrt{\lambda_{1}}} \ \ \ 0 \\ \ \ \ 0 \ \ \ \frac{1}{e}\end{bmatrix} \begin{Bmatrix}\lambda_{1} \ 0 \\ 0 \ \ 0 \end{Bmatrix} \begin{bmatrix}\frac{1}{\sqrt{\lambda_{1}}} \ \ \ 0 \\ \ \ \ 0 \ \ \ \frac{1}{e}\end{bmatrix} \begin{bmatrix}\xi\\\eta \end{bmatrix}+\left \lfloor 0 \ e \right \rfloor \begin{bmatrix}\frac{1}{\sqrt{\lambda_{1}}} \ \ \ 0 \\ \ \ \ 0 \ \ \ \frac{1}{e}\end{bmatrix} \begin{bmatrix}\xi\\\eta \end{bmatrix}+g=0$$

Simplifying the above equation further, we obtain

$$\xi^2- \eta+g=0$$ ---(7)

Step(3)- Second Translation

Now define $$\begin{bmatrix}\bar{\xi}\\ \bar{\eta}\end{bmatrix}= \begin{bmatrix}\xi \\ \eta-g \end{bmatrix}$$

This is a valid translation because $$\lambda_{1}$$ which is a function of the coefficients A,C (as defined in the equation of conics) is a constant value. Hence g which is a function of $$\lambda_{1}$$ and f is also aconstant. Thus all we are doing in the second translation is translate the $$\bar{\bar{y}}$$ by a constant. Hence making this substitution in (7) gives us the final canonical form (2)

$$\bar{\xi}^{2}- \bar{\eta} = 0$$

Lecture#38
The General form for the Equation of Conics is given below:

$$\lambda_{1}\bar{\bar{x}}^2+\lambda_{2}\bar{\bar{y}}^2+f-\frac{d^2}{4\lambda_{1}}-\frac{e^2}{4\lambda_{2}}=0$$ ---(1)

Define $$ g= -(f-\frac{d^2}{4\lambda_{1}}-\frac{e^2}{4\lambda_{2}})$$. Then (1) reduces to the form

$$\lambda_{1}\bar{\bar{x}}^2+\lambda_{2}\bar{\bar{y}}^2= g$$--(2)

Here $$\lambda_{1}$$ and $$\lambda_{2}$$ denote the eigenvalues obtained by solving the eigenvalue problem for the original equation of conics: $$Ax^2+Bxy+Cy^2+\bar{d}x+\bar{e}y+\bar{f}=0$$

Equation (2) can be further reduced to equations representing the canonical forms for an ellipse, hyperbola and parabola.

Case (i):Ellipse: $$\lambda_{1}$$ and $$\lambda_{2}$$ >0 $$\Rightarrow$$ g>0

The canonical form representing the ellipse then reduces to: $$\xi^2+\eta^2=1$$ --(3)

Case (ii): Hyperbola: $$\lambda_{1}$$ and $$\lambda_{2}$$ <0 $$\Rightarrow g\neq 0 $$

The canonical form representing the hyperbola takes two different forms depending on the sign of 'g'.

Case ii(a): g>0 $$\to \xi^2-\eta^2=+1$$ ---(4a) (This is also true for $$\lambda_{1} < 0$$,$$\lambda_{2} > 0$$ and g>0)

Case ii(b): g<0 $$\to \xi^2-\eta^2=-1$$ ---(4b) (This is also true for $$\lambda_{1} < 0$$,$$\lambda_{2} > 0$$ and g<0) Case ii(c) $$\bar{\xi}\bar{\eta}=\pm 1$$--(4c)

This is a different representation of the canonical form for hyperbolas. It is derived from either of the equations (4a) or (4b)

Note that $$[\xi,\eta]$$ represent a different coordinate system using which (2) has been transformed into the various canonical forms shown. The derivation of forms(3),(4a) and (4b) are shown in separate sections. The derivation of form (4c) from (4a)or (4b) is shown below.

Define a new coordinate system $$[\bar{\xi}, \bar{\eta}]$$ in terms of the old coordinate system $$[\xi, \eta]$$ such that $$ \begin{Bmatrix}\bar{\xi} \\ \bar{\eta}\end{Bmatrix}= \begin{Bmatrix}\xi-\eta \\ \xi+\eta\end{Bmatrix}$$

$$ \Rightarrow \begin{Bmatrix}\bar{\xi}\\ \bar{\eta}\end{Bmatrix}= \begin{bmatrix}\ \ 1 \ -1 \\ -1 \ -1\end{bmatrix}\begin{Bmatrix}\xi-\eta \\ \xi+\eta\end{Bmatrix}$$

$$\Rightarrow \begin{Bmatrix}\xi \\ \eta \end{Bmatrix}= \begin{bmatrix}\ \ 1 \ -1 \\ -1 \ -1\end{bmatrix}^{-1}\begin{Bmatrix} \bar{\xi} \\ \bar{\eta} \end{Bmatrix}= \frac {1}{2}\begin{bmatrix}\ \ 1 \ -1 \\ -1 \ -1\end{bmatrix}\begin{Bmatrix}\bar{\xi} \\ \bar{\eta} \end{Bmatrix}$$

Eqn.(4a) or (4b) can be re-written in the form: $$\begin{Bmatrix}\xi \ \eta \end{Bmatrix}\begin{bmatrix}1 \ \ \ \ 0 \\ 0 \ -1 \end{bmatrix}\begin{Bmatrix}\xi \\ \eta \end{Bmatrix}= \pm 1$$--(5)

Rewriting (5) in the new coordinate system $$[\bar{\xi}, \bar{\eta}]$$, we get:

$$ \begin{Bmatrix}\bar{\xi} \ \bar{\eta} \end{Bmatrix}\frac{1}{2}\begin{bmatrix} \ \ 1 \  -1 \\ -1 \ -1 \end{bmatrix}\begin{bmatrix}1 \ \ \ \  0 \\ 0 \ -1 \end{bmatrix}\frac{1}{2}\begin{bmatrix} \ \ 1 \   -1 \\ -1 \ -1 \end{bmatrix}\begin{Bmatrix}\bar{\xi} \\ \bar{\eta} \end{Bmatrix}= \pm 1$$

$$ \Rightarrow \frac{1}{4} \begin{Bmatrix}\bar{\xi} \ \bar{\eta} \end{Bmatrix}\begin{bmatrix} \ \ 1 \  -1 \\ -1 \ -1 \end{bmatrix}\begin{bmatrix}1 \ \ \ \  0 \\ 0 \ -1 \end{bmatrix}\begin{bmatrix} \ \ 1 \   -1 \\ -1 \ -1 \end{bmatrix}\begin{Bmatrix}\bar{\xi} \\ \bar{\eta} \end{Bmatrix}= \pm 1$$

$$\Rightarrow \frac{1}{4}\begin{Bmatrix}\bar{\xi} \ \bar{\eta} \end{Bmatrix}\begin{bmatrix}\ \ \ 0 \ -2 \\ -2 \ \ \ \ 0 \end{bmatrix}\begin{Bmatrix}\bar{\xi} \\ \bar{\eta} \end{Bmatrix}= \pm 1$$

$$\Rightarrow \bar{\xi}\bar{\eta}=\pm 1$$ -(4c)

Thus we see that by a rotation of coordinates, form (4a) or (4b), which is the standard canonical form representign a hyperbola reduces to form (4c). This simplified form is used mostly in algebraic calculations and even in high-school geometry lessons.

Term Project: Understanding the Plasma-Wall Problem
I am a doctoral candidate enrolled in the Mechanical & Aerospace Engineering Department at the University of Florida, Gainesville. I work in the Computational Plasma Dynamics Laboratory & Testing Facility under the supervision of Dr.Subrata Roy.



On earth, we know and have learnt to live with only three states of matter. However Sir William Crookes identified Plasma, the fourth state in 1879. It is by far the most common form of matter. It is present in the stars and in the tenuous space between them. It is present, though invisible, in almost 99.9% of the visible universe. Ideally plasma consists of neutrals, electrons and ions. However commercial plasma is ionized gas. In order to separate electrons and ions from atoms, energy is needed and this energy can be in any form, thermal, electric etc. Most commonly, a gas is passed through a pair of electrodes and an electric field. This electric field provides the energy needed to ionize the gas into atoms and ions. Plasma has helped create a better understanding of the world we live in. It provides us the prospect of abundant energy. It provides alternatives to surface cleaning, waste removal, in the medical industry and many more.

In my research, we focus on Dielectric Barrier Discharges(DBDs) & their applications. Commerically produced plasma is classified into two types- Volume Plasma (wherein the plasma is produced in the discharge space between two electrodes) & Surface Plasma (wherein both electrodes are embedded on the same surface). DBDs are a kind of surface plasma. They are generated between two electrodes with a dielectric barrier in between them. The electrodes are separated by a gap of a few millimetres. Typically DBDs are operated at voltages of 1-100 kV and frequencies of 50Hz-1MHz. When electric field is sufficiently high to cause breakdown of the discharge gas, a large number of micro-discharges can be observed emanating from the electrodes.

Current Research Focus- Plasma Sterilization Plasma Sterilization is a process in which low-temperature plasma is utilized to kill micro-organisms. Our research focus is on using Dielectric Barrier Discharges for this purpose. Plasma is produced usually by passing a gas through an electric field and ionizing it. Dielectric Barrier Discharges are a special case, in which two electrodes are embedded into the top and bottom of a dielectric surface (an insulator). An electric field is produced between these two electrodes. When a gas is passed through this electric field, it is ionized, producing ions, electrons, neutrals and UV photons. All of these species collide with the micro-organism, reacting with its essential elements and forming harmless compounds. The micro-organism, when devoid of its essential elements becomes harmless and is thus inactivated. More information about preliminary experiments on this subject is given here.

Abstract
One of the oldest problems that has sparked a lot of research in plasma physics is the plasma-wall problem. The sheath is a region of electric field that exists between the plasma and the boundary into which it runs. The basic problem of plasma running into a wall needs to be understood in order to better visualize some of the other aspects of plasma physics. In this paper, we attempt to understand the processes happening in the sheath region. There exist a set of governing equations, akin to the governing equations of fluid dynamics that outline the physical processes running in the plasma. These equations are discretized using a simple, first order finite difference scheme and then solved using the Newton-Raphson Method, as described in a previous paper[3]. The boundary conditions applied however are those pertaining to the sheath region existing between plasma and a solid wall. The objective in conducting this analysis is to investigate the accuracy of the numerical scheme used. If the numerical scheme is accurate, the electron density at the wall should go to zero. Further a probability density function is constructed using electron density data generated from the solver code.

Introduction
To completely describe the state of plasma, we would need to write down all the particle locations and velocities and describe the electromagnetic field in the plasma. But considering the zillions of particles, this would be an impossible task! In order to make our life easier, we construct models which simplify the processes taking place in plasma.

Models are of two types:
 * Fluid Model- This describes plasmas in terms of smoothed quantities like density and averaged velocity around each position. A general description is the two-fluid picture, where ions and electrons are described separately. Fluid models usually model the particles using a Maxwell-Boltzmann Distribution. However they are usually unable to capture accurately finer structures such as beams or double layers.
 * Kinetic Model- While being more accurate than the fluid model, it is usually more computationally intensive. It describes the particle velocity distribution function at each point in the plasma and therefore does not need to assume a distribution function like the fluid model. The more commonly used Particle-In-Cell (PIC) technique includes kinetic information by following the trajectories of a large number of individual particles

Understanding plasma sheaths is perhaps one of the oldest problems in plasma physics. A plasma sheath is the localized electric field that separates plasma from a material boundary. It confines the more mobile species (electrons) in the plasma and accelerates the less mobile species (ions) out of the plasma and toward the walls. This means that in the typical case, where the electrons are more mobile than the positively charged ions, the electric field in the sheath points towards the boundary.This behaviour can be further influenced by external processes such as variations in the applied electric field or the imposition of an external magnetic field.

Modeling of the behaviour of species in the Sheath Region
Many numerical models have been developed to describe Sheaths. The simplest visualization of this process is the behaviour of electrons and ions in the sheath region. Consider a planar slab of collisionless plasma, modeled in 1D space, using a fluid approach and assuming a Boltzmannian electron distribution. The governing equations for such a model are given below:

Continuity Equation for Ion flux:

$$ \frac{\partial (n_{i}v_{i})}{\partial x}= zn_{e}$$ ---(1)

Continuity Equation for Electron flux:

$$ \frac{\partial (n_{e}v_{e})}{\partial x}= zn_{e}$$ ---(2)

Momentum of ions:

$$ \frac{\partial (Mn_{i}v_{i}^{2})}{\partial x}+en_{i}\frac{\partial(\psi)}{\partial{x}}= 0 $$ ---(3)

Momentum of electrons:

$$ \frac{\partial (mn_{e}v_{e}^{2}+kn_{e}T_{e})}{\partial x}-en_{e}\frac{\partial(\psi)}{\partial{x}}= 0$$(4)

Poisson Equation:

$$ \epsilon_{0}\frac{\partial^2 \psi}{\partial x^2}= -e(n_{i}-n_{e})$$ --(5)

where k is the Boltzmann Constant, m is the electron mass, M is the ion mass, Te is the electron temperature,ve is the electron velocity, viis the ion velocity, $$\psi$$ is the electric potential, z is the frequency of direct ionization, e is the electron charge, $$\epsilon_{0}$$ is the permittivity of free space, niis the ion density and ne is the electron density.

More details on these equations and their non-dimensionalization can be obtained from the paper by Sternberg and Godyak on Asymptotic Matching and the Sheath Edge

For our purpose, we take the normalized versions of Equations (1)-(5) used in the paper by Roy & Gaitonde.

$$\frac{\partial{(y_{i}V_{i}})}{\partial{x}}= y_{e}= exp(-\phi)$$--(5)

$$ V_{i}\frac{\partial{V_{i}}}{\partial{x}}-E+V_{i}\frac {y_{e}}{y_{i}}=0$$---(6)

$$ \frac{\partial{E}}{\partial{x}}- \kappa^{2}(y_{i}-y_{e})= 0$$ --(7)

$$ E= \frac{\partial{\phi}}{\partial{x}}$$(8)

where $$y_{i}$$ is the normalized ion density, $$y_{e}$$ is the normalized electron density, $$V_{i}$$ is the normalized ion velocity, E is the electric field strength and $$\phi$$ is the electric potential.

Solution through Numerical Schemes
A finite difference scheme is applied to discretize the first derivatives in Equations(5)-(7). Equation(8) is just needed to calculate the electric potential $$\phi$$ Finite difference schemes are usually applied for complex, non-linear partial differential equations, wherein the solution space is divided into n number of grid points and the required value at each grid point is calculated. For the single derivative, there are typically three types of schemes: Forward difference, Backward difference and Central Difference. Each of these have a different order of accuracy and are applied depending on how accurate the numerical scheme applied has to be.

Forward Difference$$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n+1}-A_{n}}{\delta x}$$

Backward Difference$$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n}-A_{n-1}}{\delta x}$$

Central Difference $$\Rightarrow \frac{\partial {A}}{\partial x}= \frac {A_{n+1}-A_{n-1}}{2\delta x}$$

where An+1 is the value of A(the quantity being calculated) at the n+1th grid point and An-1 is the value of A at the n-1th grid point.

The Central Difference Scheme usually has an order of accuracy higher than the other two schemes but this does not necessarily mean that using this scheme might give us the most accurate solution. The nature of the PDE being solved should also be taken into account. This is where we define terms like 'stability' and 'consistency' of a finite difference scheme. For more information, Computational Fluid Dynamics- John.D.Anderson is a very informative book to read. The basics of fundamental numerical schemes as well as more complex ones are outlined along with illustrative applications to common problems in Fluid Mechanics like the Burger's Equation.
 * Hence once Equations (5)-(7) are discretized using forward difference, they reduce to a set of non-linear equations in which the unknowns are yi, Vi and E. This set of non-linear equations is solved using the Newton-Raphson Method. Once these are known, ye can be calculated using Equation(5).

Newton Raphson Method

Consider a typical problem that gives N functional relations to be zeroed involving variables xi, i=1,2....N.

f(xi)= 0, i=1,2....N

Each of these functions can be expanded into a Taylor Series i.e $$ f(x_{i}+\delta x_{i})= f(x_{i})+ \delta x_{j}\frac {\partial{f_{i}}}{\partial{x_{j}}}+O(\delta x^2)$$

Assuming $$f \to 0$$, we get $$ \Sigma \alpha_{ij}\delta x_{ij}= \beta_{ij} $$ where $$ \beta_{ij}= -f_{i} $$

$$\alpha_{ij} $$ is known as the Jacobian. Here Equations 5-7 are the equations that are discretized and contain three unknowns. So $$\alpha$$ is going to be a 3x3 matrix and $$\beta$$ will be a 3x1 matrix. $$\delta x_{ij}$$ is then solved for iteratively in MATLAB.


 * Now, once the electron density distributions and ion-density distributions are calculated, how do we estimate what happens down the line, say after like 100 iterations? Does the electron density in the Sheath Region remain zero always, as it should? Does the ion-density increase or decrease? This check achieves two purposes: The accuracy of the numerical scheme is verified and the sheath behaviour after a long time is understood. This check is performed with the help of a probability density analysis. A MATLAB code was written to simulate the entire process described above. This code furthe takes the electron-density data generated, divides it into 100 blocks and constructs a probability density function. This is repeated for every iteration so that we know the nature of the pdf after about 100 iterations.

Results


Fig(1) shows the published results from the paper by Roy et.al[3]. In this paper, equations (1-7) to (1-10) are solved using a second-order time accurate, implicit finite element method using a computational domain extending from x=0 to the wall. Fig. 1(a) is a plot of ion density versus position x. Fig. 1(b) is the plot of ion velocity as a function of position, x. Fig. 1(c) is the plot of electric potential as a function of x, while Fig. 1(d) represents the plot of electric field versus x. As expected, since the more mobile species is confined within the plasma and the less mobile species is accelerated out, it is seen accordingly that at xw~0.71, the ion velocity increases. Since electrons are being confined, one would expect electric field at wall to rise (more charge being accumulated). This is seen in Fig. 1(c) and consequently in Fig. 1(d).

In Fig.(2), the plasma parameter distributions obtained from the MATLAB code used are given. In fig. 2(a), as can be seen, ion density prediction matches the trend as predicted in Fig.1(a) i,e ion density decreases in the direction of the wall. This is expected. Fig. 2(b), which shows the plot of ion velocity versus x also seems to match the trend predicted in 1(b). It makes senses too because one would expect ion velocity to increases as sheath approaches. Fig.2(c) makes perfect sense because electron density decreases to zero as wall is approached. The value of x at which ne goes to zero is a little earlier than expected (at x=0.2). Also to be noted is the important fact that ion density at wall does not go to zero while electron density does go to zero, which simply confirms the fact the electrons stay behind in the plasma while ions accelerate towards the wall. Fig. 2(d) is a plot of the potential versus x, which also matches the plot shown in fig.1(c). This tells us that the code is predicting the behaviour of the species in the sheath correctly to some extent. The linear trend present in fig. 2(a), 2(b), 2(d) is of slight concern as we are looking for more of an exponential trend. This could be because of an error in the code or due to slight errors introduced by the numerical scheme used. Note that the results in fig.(1) are the result of a finite element code while the results in fig. (2) are those computed from a finite difference scheme. Issues of accuracy tend to exist.

Fig. 3 shows the pdf analysis for the electron density. Electron density data is taken from the code, divided into 100 blocks and a pdf is constructed. Fig.3(a) shows this pdf as a bar diagram. Further the convergence of this pdf is investigated by calculating moments as high as Ethumb|left|Fig.3:PDF analysis for electron density dataNavya4.jpg There are a couple of important points to note in these two figures. Fig.3(a) shows that there is a 4% probability of electron density going to zero at the wall. Fig. 4(a) shows that the ion density does not go to zero at the wall. Fig. 3(b) and Fig. 4(b) confirm convergence of the pdf for both electron density and ion density data. Note also that the electron density is steady after about 40 iterations while the ion density is steady after about 20 iterations

Conclusion & Future Work
The governing equations for the plasma wall problem were solved using a finite difference scheme to discretize them. The electron and ion density are calculated using a MATLAB code. The computed distributions of ion and electron density exhibit the trend shown in the published results as shown in Fig.(1). The linear trend may be due to accuracy errors in solver code. However the ion and electron density plots are very similar to ion and electron behaviour in sheath region. Also the PDF analysis shown in Fig. (3) and (4) for electron and ion density data respectively shows that there is convergence after 100 iterations, proving that the results predicted by the code are accurate. Also we see that there is a 4% probability of the electron density going to zero. This value might further reduce, provided the accuracy errors in the code are resolved.

Future work involves using more complex, higher order numerical schemes to discretize the governing equations. Also a simplistic case of the plasma-wall phenomenon was examined here. More complex phenomenon can be introduced into this simulation in order to further understand this phenomenon.