1

I'm trying to solve a system of two equations, both of which use the same interpolating function. If the background would help, I'm trying to find the depletion width of a p-n junction after I've evaporated a dopant (donor) onto it and I've let it diffuse for a while. I'm trying a simple case first, where the semiconductor has a constant acceptor concentration of .8 atoms/cm^-3 and the donor concentration is some interpolation function that I got by solving the diffusion equation using the Crank-Nicolson method. I find the junction (where the two concentrations are equal) and then I solve the two simultaneous equations: 1.) Derivative of potential is zero at both boundaries 2.) Continuity of the scalar potential @ the junction.

Here's my code (just signed up for Mathematica SE, please bare with me):

Solve[{0.8*lp == (solution[ln, 5])*ln, (-solution[0, 5])*(ln)^2 == 0.8*(lp)^2}, {lp, ln}]

where 'solution' is the interpolation function (the donor concentration in the semicondcutor). However, when I try and solve this system of equations it gives me the error: "Solve::inex: Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help."

I'm not sure if this is something that Mathematica can even do, but any help would be greatly appreciated.

Edit: Here's my code for how I got the interpolation function, 'solution

M = 100; Δt = 0.05; Δx = 0.1; junc = 3.16294;
ρ0[x_] := N[Erfc[x]]; 
ρ[0, n_] := ρ[1, n]; 
ρ[M, n_] := ρ[M - 1, n]; 
sol[0] = Table[ρ[j, 0] -> ρ0[j*Δx], {j, 1, M - 1}];
sol[n_] := sol[n] = Module[{vars, eqns}, vars = Table[ρ[j, n], {j, 1, M     - 1}]; 
 eqns = Table[ρ[j, n] - ρ[j, n - 1] == (Δt*(ρ[j + 1, n] - 2*ρ[j, n] + ρ[j - 1, n] + ρ[j + 1, n - 1] - 2*ρ[j, n - 1] + ρ[j - 1, n - 1]))/(2*Δx^2), 
    {j, 1, M - 1}] /. sol[n - 1]; Solve[eqns, vars][[1]]]
solution = Interpolation[Flatten[Table[Table[{Δx*j - junc, n*Δt, ρ[j, n]}, {j, 0, M}] /. sol[n], {n, 0, 100}], 1]]
george2079
  • 38,913
  • 1
  • 43
  • 110
hijasonno
  • 41
  • 2
  • Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! –  Nov 16 '16 at 06:55
  • Have a look at How to get rid of warnings when using Solve on an equation with inexact coefficients? –  Nov 16 '16 at 06:59
  • 1
    We can't run your code, because you haven't provided solution and junc. – Feyre Nov 16 '16 at 10:31
  • Just edited with code for 'solution' interpolation function. – hijasonno Nov 16 '16 at 18:21
  • A quick look at solution shows that $f\left(ln,lp\right) \ge 0$ in the domain of the data used in its construction. That means that the first equation will only be solved for $\left(ln,lp\right) = \left(0,0\right)$, which is not a solution to the second equation. – Marchi Nov 16 '16 at 18:44
  • Oops, I don't think that negative sign should be there in the first set of equations, it should be there in the second set however. – hijasonno Nov 16 '16 at 18:53
  • So I think what I actually have to do is shift the x-axis so that 'junc' sits at x = 0 in order for this to give me any accurate information. The way I've tried to solve Poisson's equation, I unknowingly implied that 'junc' is equal to zero which I believe is causing all of these problems. – hijasonno Nov 16 '16 at 19:01
  • I've fixed the code to reflect the shifting in origin, but I still get the same error. – hijasonno Nov 16 '16 at 19:49
  • 2
    Try using FindRoot. Using the most updated information I'm aware of the solution to the system is $\left(ln,lp\right) = \left(0,0\right)$. Try FindRoot[{0.8*lp == (solution[ln, 5])*ln, (-solution[0, 5])*(ln)^2 == 0.8*(lp)^2}, {{lp, 1}, {ln, 1}}] – Marchi Nov 16 '16 at 20:34
  • solution turns out to be positive everywhere so you can readily see the second equation of the form -solution (ln)^2 == (lp)^2 can only have the trivial 0,0 solution. – george2079 Nov 16 '16 at 20:52
  • @Marchi, thanks a lot, that worked. I guess the trick is for interpolation functions is to use built-in functions that evaluate the interpolation function at each step rather than using something like 'Solve,' which doesn't do that. – hijasonno Nov 17 '16 at 05:35

0 Answers0