Iterative Solution of Highly Nonlinear Differential Equations Using 
Mathematica

Frank Kampas

Lockheed Martin*

The  work presented here is not related to the author's  work for Lockheed 
Martin

Introduction

Nonlinear boundary value differential equations are usually solved with  the "shooting method".  In this technique, the initial conditions are  adjusted until the boundary conditions at the other boundary are satisified.   In situations in which the shooting method fails, the iterative  "relaxation" method can be used.  Initial guesses at the solution are  improved repeatedly.  The relaxation method can be implement using  "quasi-linearization".  The differential equation is linearized about the  guessed solution and the linearized boundary value problem is solved.  This  technique can be implemented in a very straight-forward fashion in Mathematica because the numerical solutions to differential equations are  interpolating functions.

Quasi-linearization is first demonstrated with a  single differential equation and then with the 5 coupled differential  equations which describe a p-n junction.

Troesch's Equation

The shooting method fails for Troesch's equation, y'' = 5 Sinh(5 y), y(0)=0, 
y(1) =1, where the solution has a singularity at approximately at x = 1.03, 
unless the initial estimate of the value of the derivative at x = 0 is very 
good.  In the following calculation, shootsoln[d_] is the solution to 
Troesch's equation, as a function of d, which is equal to y'(0).  FindRoot is 
used to find the value of d for which the solution satisfies the boundry 
condition at x = 1.

shootsoln[
        d_?NumberQ] := NDSolve[{y''[x] - 5 Sinh[5 y[x]] == 0, 
            y[0] == 0, y'[0] == d}, 
          y, {x, 0, 
            1}][[1]]

yd[d_?NumberQ] := y[1] /.  shootsoln[d]

\(soln = 
      y\  /. \ \((shootsoln[
              d]\  /. 
		\(FindRoot[
                yd[d] == 1, {d, 0, 
                  1}]\)[[1]])\)\)

In the previous example, y'(0) was estimated as being between 0 and 1.0.  
However, if it were known that y'(0) is between -0.05 and 0.05, then the 
shooting method can be employed.

\(soln = 
      y\  /. \ \((shootsoln[
              d]\  /. 
		\(FindRoot[
                yd[d] == 1, {d, \(-0.05\), 
                  0.05}]\)[[1]])\)\)

Note the correct value of y'[0].

soln'[0]

When the shooting method fails, quasi-linearization can be employed.  The 
non-linear differential equation is linearized about a function and the 
linearized differential equation is solved using the shooting method.  This 
is repeated until the solution has sufficient accuracy.

The following function linearizes the left-hand side of the differential 
equation by implementing a Taylor series to first order.  The functions and 
the derivatives in the expression are input as vars, and the base functions 
and derivatives about which expr is linearized are input as noms.

lin[expr_, vars_, noms_] := 
      Module[{r = 
              Thread[vars -> noms]}, (expr  /.  
                r) + (vars - 
                  noms) . ((D[expr, #] &  /@  vars)  /.  
                  r)] // Simplify

Troesch = y''[x] - 5  Sinh[5  y[x]];

linTroesch[ys_] = 
      lin[Troesch, {y[x], y''[x]}, {ys[x], ys''[
            x]}]

The shooting solution is now a function of the starting function ys and of 
y'(0).

Remove[shootsoln,  yd]

shootsoln[ys_, 
        d_?NumberQ] := NDSolve[{linTroesch[ys] == 0, 
            y[0] == 0, y'[0] == d}, 
          y, {x, 0, 
            1}][[1]]

yd[ys_, d_?NumberQ] := 
      y[1] /.  shootsoln[ys, d]

The function iter improves the solution using the shooting method to solve 
the linearized differential equation.  Note that its input and output are 
anonymous functions.  This allows the evaluation of derivatives.

iter[ys_] := 
      y  /.  ( 
          shootsoln[ys, 
              d]  /.  FindRoot[
                yd[ys, d] == 
                  1, {d, ys'[
                    0],  ys'[
                      0] +  001}][[1]])

The function #&, which satisfies the initial conditions, is used as the 
starting function.

l = NestList[iter, # &, 7];

The following plot shows the convergence

Plot[Evaluate[Through[l[x]]], {x, 0, 1}];

All of the solutions satisfy the boundary conditions.

Through[l[0]]

Through[l[1]]

The value of y'(0) converges to the value found earlier by the shooting 
method

Table[l[[i]]'[0], {i, 8}]

P-N Junction Transport Equations

The equations describing the electron and hole transport in a p-n junction 
cannot be solved by the shooting method, but can be solved with 
quasi-linearization.

The two following equations give the electron density n as a function of  potential v and the electron quasi-Fermi level ϕn and the hole density p  as a function of the potential and the hole quasi-Fermi level ϕp.  ni is  the intrinsic carrier density, β = (1/(k T)) where k is Boltzmann's constant and T is the absolute temperature.

n[v_, ϕn_] = ni*Exp[β (v - ϕn)];
p[v_, ϕp_] = ni*Exp[β (ϕp - v)];

The total charge density chgDens[x] is equal to the dopant charge density 
dp[x] + the electron and hole charge densities.

chgDens[x_] = 
        dp[x] - n[v[x], ϕn[x]] + p[v[x], ϕp[x]];

The following 5 terms are equal to 0 when the differential equations are 
satisfied.  jn and jp are the electron and hole current densities, 
respectively.  The first two equations give the relationship between the 
current densities and the quasi-Fermi levels.  The third equation is 
Poisson's equation.  The fourth and fifth equations reflect the fact that the 
current densities are independent of x in steady-state.

eq1 = ϕn'[x] + 
          jn[x]/(q μn n[v[x], ϕn[x]]);
eq2 = ϕp'[x] + 
          jp[x]/(q μp p[v[x], ϕp[x]]);
eq3 =  v''[x] + (q/ϵ) chgDens[x];
eq4 = jn'[x];
eq5 = jp'[x];

eqs is a list of the five equation left-hand sides.

eqs = {eq1, eq2, eq3, eq4, eq5};

The linearized differential equations are functions of the starting 
functions, whose names end in s.

delin[vs_, ϕns_, ϕps_, jns_, jps_] = 
      lin[eqs, {v[x], v'[x], v''[x], ϕn[x], ϕn'[
              x], ϕp[x], ϕp'[x], jn[x], 
            jp[x]}, {vs[x], vs'[x], vs''[x], ϕns[
              x], ϕns'[x], ϕps[x], ϕps'[x], jns[x], 
            jps[x]}] // FullSimplify

The following are the initial conditions for the linearized equations.  e0  is the electric field as x = 0.  The quasi-Fermi levels are equal to each  other and to ϕ0 at the left boundary.

ics[e0_, jn0_, 
          jp0_] = {v[0] == 
            v0, ϕn[0] == ϕ0, ϕp[0] == ϕ0, -v'[0] == 
            e0, jn[0] == jn0, jp[0] == jp0};

The following are values of various constants

\(\(\(q = 1.6\ 10\^\(-19\)\ ; \  (*\ 
      the\ charge\ on\ an\ electron\ in\ Coulombs\ *) 

    k = 1.38\ 10\^\(-23\); \  (*\ Boltzmann' s\ constant\ *) 

    t = 300\ ; \  (*\ the\ temperature\ in\ degrees\ Kelvin\ *) 
β = 
      q\/\(k\ t\);

    \(\(μn = 1500\) \); \  (*\ 
      the\ electron\ mobility\ *) 
\(\(μp = 450\) \); \  (*\ 
      the\ hole\ mobility\ *) 
\(\(ni = 1.45\ 10\^10\) \); \  (*\ 
      the\ intrinsic\ carrier\ density\ *) 
\(\(ϵ = 
        11.9*8.85\ 10\^\(-14\)\) \); \  (*\ 
      the\ permitivity\ of\ silicon\ *) 
\(\(l = 
        10\^\(-4\)\) \); \)\(\ \)\( (*\ 
      the\ total\ thickness\ of\ the\ junction\ *) \)\)\)

The charge density from the dopants is -dop for 0 < x < (l/(2 ))and dop for (l/2) < x <  x.  Therefore the junction is a p-n junction since the negatively charged  acceptors are between 0 and (l/2) and the positively charged donors are between (l/2) and l.

dp[x_] = dop*Sign[x - (l/2)];

The following are the electric potential across the junction and the dopant 
concentration.

\(\(\(\(dop = 10\^15\) \)\(;\)\(\ \)\( (*\ 
      the\ dopant\ concentration\ *) \)\(
\)\(\(v0 = 
        0.0\) \)\(;\)\(\ \)\( (*\ 
      the\ potential\ at\ the\ left\ boundary\ *) \)\(
\)\(\(v1 = 
        0.45\) \)\(;\)\(\ \)\( (*\ 
      the\ potential\ at\ the\ right\ boundary\ *) \)\)\)

The built-in potential can be calculated from the concentration of dopants.  
The fact that the potential difference across the junction is less than the 
built-in potential indicates that the junction is forward biased.

vbi = (β^(-1)) Log[((dop^2)/(ni^2))]

At the boundaries, the quasi-Fermi levels are equal and the total charge 
density dp - n + p is 0.  This enables calculation of the values of the 
quasi-Fermi levels at the boundaries.

ϕ0 = ϕ0  /.  FindRoot[
            dp[0] - n[v0, ϕ0] + p[v0, ϕ0] == 
              0, {ϕ0,  3}][[1]]

ϕ1  =  ϕ1  /.  FindRoot[
            dp[1] - n[v1, ϕ1] + p[v1, ϕ1] == 
              0, {ϕ1,  3}][[1]]

The following function linsoln gives the solution to the linearized  equation as a function of the starting functions vs, ϕns, ϕps, jns  and jps, and as a function of the initial conditions e0, jn0 and jp0.

\(linsoln[vs_, ϕns_, ϕps_, jns_, jps_, e0_?NumberQ, 
        jn0_?NumberQ, 
        jp0_?NumberQ] := 
	\(NDSolve[
		Thread[
              delin[vs, ϕns, ϕps, jns, jps] == 0]~Join~
            ics[e0, jn0, jp0], {v, ϕn, ϕp, jn, jp}, {x, 0, l}, 
          Compiled -> False, 
          InterpolationPrecision -> 
            5]\)[[1]]\)

bcvals gives the values at the right-hand boundary which need to be satisfied.  It is similar to the function yd in the earlier example.

bcvals[vs_, ϕns_, ϕps_, jns_, jps_, e0_?NumberQ, 
        jn0_?NumberQ, jp0_?NumberQ] := {v[l], ϕn[l], ϕp[l]}  /. 
        linsoln[vs, ϕns, ϕps, jns, jps, e0, jn0, 
          jp0]

The function iter calculates improved functions for the potential, the 
electron and hole quasi-Fermi levels, and the electron and hole current 
densities, using linsoln.  The starting values for the required initial 
conditions are calculated from the functions.

\(iter[{vs_, ϕns_, ϕps_, jns_, jps_}] := {v, ϕn, ϕp, 
          jn, jp}\  /. \((linsoln[vs, ϕns, ϕps, jns, jps, e0, jn0, 
              jp0] /. \ FindRoot[

              bcvals[vs, ϕns, ϕps, jns, jps, e0, jn0, 
                  jp0] == 
				{v1, ϕ1, ϕ1}, {e0, \(-\  .95\)\ \(vs'\)[0], \(-1.05\)\ \(vs'\)[0]}, {jn0,  .95\ \ jns[0], 
                1.05  jns[0]}, {jp0,  .95\ jps[0], 1.05\ jps[0]}, 
              MaxIterations -> 6])\)\)

The initial potential and quasi-Fermi level functions are assumed to be 
linear.  The initial current density functions are assumed to be constant.

i0 = {(v0 + ((v1 - 
                        v0)/l) #) &, (ϕ0 + ((ϕ1 - ϕ0)/l) #)  &, (ϕ0 + ((ϕ1 - ϕ0)/l) #)  &, 
          (10^(-7)) &, (10^(-7))  &};

The function iter is applied 7 times.  The number of MaxIterations for 
FindRoot was limited because it is the repeated application of iter which 
improves the solution.  Therefore precision warning messages are generated.

iterations = NestList[iter, i0, 7];

The following plot shows the convergence of  v[x] .

Plot[
        Evaluate[Table[iterations[[nn, 1]][x], {nn, 7}]], {x, 0, 
          l}];

The next plot shows the convergence of ϕn[x]  and ϕp[x] .

Plot[
        Evaluate[
          Join[Table[iterations[[nn, 2]][x], {nn, 7}], 
            Table[iterations[[nn, 3]][x], {nn, 7}]]], {x, 0, l}, 
        PlotRange -> {{0, l}, { .16,  .30}}];

This table shows the convergence of the the electron current density jn

Table[iterations[[nn, 4]][0], {nn, 7}]

The convergence of jp.

Table[iterations[[nn, 5]][0], {nn, 7}]

The convergence of the electric field at 0, e0 .

-Table[iterations[[nn, 1]]'[0], {nn, 7}]

The converged potential, electron quasi-Fermi level, hole quasi-Fermi  level, electron current, and hole current are given by vf, ϕnf,  ϕpf, jnf and jpf.

{vf, ϕnf, ϕpf, jnf, jpf} = 
        iterations[[7]];

The total charge density can be calculated from the second derivative of the 
potential, per Poisson's equation.  Most of the junction is depleted of 
carriers.

Plot[((-ϵ)/q) vf''[x], {x, 0, l}, 
        PlotRange -> {-(10^15), (10^15)}];

The electron and holes densities can be calculated from the potential and the 
quasi-Fermi levels.

Plot[{n[vf[x], ϕnf[x]], p[vf[x], ϕpf[x]]}, {x, 0, l}, 
        PlotRange -> All];

Conclusion

Complex coupled nonlinear boundary value differential equations can be  solved in Mathematica using quasi-linearization by taking advantage of Mathematica's unique use of interpolating functions for the numerical solution of  differential equations.