2
$\begingroup$

I'm trying to solve the drift-diffusion equation $$\frac{\partial\rho}{\partial t} = \frac{\partial^2 \rho}{\partial x^2} + \frac{\partial \rho}{\partial x} \rho(x,t)$$ using a finite difference method, where I solve the diffusion equation for half of the spatial grid, and this equation for the second half of the spatial grid. Here's my code:

M = 100; Δt = 0.05; Δx = 0.1; 
ρ0[x_] := N[Erfc[x]]; 
ρ[0, n_] := ρ[1, n]; 
ρ[M, n_] := ρ[M - 1, n]; 

sol[0] = Join[Table[ρ[j, 0] -> ρ0[j*Δx], {j, 1, 50}], 
Table[ρ[j, 0] -> ρ0[j*Δx], {j, 51, M - 1}]];

sol[n_] := sol[n] = Module[{vars, eqns}, vars = Table[ρ[j, n], {j, 1, M - 1}]; 
 eqns = Join[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, 50}], 
    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) + 
       (Δt/Δx)*(j*Δx*(ρ[j, n] - ρ[j, n-1])ρ[j, n]), {j, 51, 99}]] /. sol[n - 1]; 
 FindRoot[eqns, ({#, 0.4} &) /@ vars]]

I'm currently applying Neumann conditions, but I want to employ Robin conditions so that $\rho$ is conserved throughout the grid. How do I go about doing that for a finite difference method?

$\endgroup$

closed as off-topic by user21, m_goldberg, corey979, MarcoB, happy fish Nov 22 '16 at 14:22

  • The question does not concern the technical computing software Mathematica by Wolfram Research. Please see the help center to find out about the topics that can be asked here.
If this question can be reworded to fit the rules in the help center, please edit the question.

  • $\begingroup$ I wonder if your question is about Mathematica the software, and not about the underlying mathematics of your problem? $\endgroup$ – MarcoB Nov 22 '16 at 4:58
  • $\begingroup$ @MarcoB, perhaps? I just don't really understand how to apply Robin boundary conditions in a finite difference setting in general. It's very clear for Neumann conditions (just take the forward and backward derivatives) and Dirichlet conditions, but I'm just lost on how to do it for Robin conditions. I was thinking that instead of setting the boundary conditions at the top (in lines 3 & 4 of the code) I could amend sol[n_] to make the eqns include the Robin conditions on the boundary, but wasn't sure if that made any sense. $\endgroup$ – hijasonno Nov 22 '16 at 6:41
  • $\begingroup$ Once again, are you practicing coding FDM, or you just want to solve the PDE no matter what method is used? $\endgroup$ – xzczd Nov 22 '16 at 7:34
  • $\begingroup$ I'm voting to close this question as off-topic because I think I should go to scicomp.stackexchange.com $\endgroup$ – user21 Nov 22 '16 at 9:56
  • 3
    $\begingroup$ I'm voting to close this question as off-topic because the issue it raises is not a Mathematica issue but a mathematics issue. That it is formulated in terms of Mathematica is not sufficient to make it an appropriate question for Mathematica.SE. $\endgroup$ – m_goldberg Nov 22 '16 at 10:33