Solving coupled PDEs numerically on a semi-infinite domain with no-flux boundary conditions

Question

I have the following system of PDEs for which I have given parameters $\gamma, \tau$ and $\mu$,

$$\begin{align} T_t = &\ \gamma\,(L +\tau F-T)\\ F_t = & -F_x-(F-LT)\\ L_t = &\ \mu L_{xx}+(F-LT)\end{align}$$

with no-flux boundary conditions at $x=0$ and $x=\infty$ for $F$ and $L$. The initial conditions are $L=T=F=e^{-x^2}$.

Mathematica failed to produce results because the order of the derivative in the boundary condition for $F$, $(F_x(0)=F_x(\infty)=0)$ is the same as the order of the derivative of $F$ w.r.t $x$ in the PDE.

Any alternatives to Mathematica?

Answer 1

Here is something I threw together, it is definitely incorrect, as the program is throwing errors. In particular, a SolutionVariableNumberError is raised. Likely this is because the equations are written incorrectly. The likely suspect, I'm guessing is in the SourceTerm's. This is a start though for how one might go about solving these equations in FiPy.

Perhaps someone in the forum has more experience with this tool and can elaborate on what is going wrong here. Other than that if you want to continue with this tool, I suggest reading and studying more of FiPy's examples. You could also post to the their mailing list, and ask for some help writing the equations. I've seen others do this on the list, as I have been following it for a while. The project maintainers are quite helpful. (I have no affiliation to the project, I'm just a fan).

from fipy import CellVariable
from fipy import DiffusionTerm
from fipy import Grid1D
from fipy import ImplicitSourceTerm
from fipy import numerix
from fipy import PowerLawConvectionTerm
from fipy import TransientTerm

mu = 1.
gam = 1.
tau = 1.

m = Grid1D(nx=100,Lx=1.)

L = CellVariable(mesh=m, hasOld=True)
F = CellVariable(mesh=m, hasOld=True)
T = CellVariable(mesh=m, hasOld=True)

x = m.cellCenters

L.value = numerix.exp(-x**2)
F.value = numerix.exp(-x**2)
T.value = numerix.exp(-x**2)

eqn0 = TransientTerm(var=L) == DiffusionTerm(mu, var=L) + (ImplicitSourceTerm(var=F) - ImplicitSourceTerm(var=L*T))
eqn1 = TransientTerm(var=F) == PowerLawConvectionTerm(-1., F) - (ImplicitSourceTerm(var=F) - ImplicitSourceTerm(var=L*T))
eqn2 = TransientTerm(var=T) == gam * (ImplicitSourceTerm(1., var=L) + ImplicitSourceTerm(tau, var=F) + ImplicitSourceTerm(-1.,var=T))

eqn = eqn0 & eqn1 & eqn2

for t in range(100):
    L.updateOld()
    F.updateOld()
    T.updateOld()
    eqn.solve(dt=1.e-3)

Answer 2

Seeing what are your boundary conditions, you should use an spectral(Galerkin) method. I am trying to solve sometrhing similar, but in 9 dimentions, And the way to treat the field fading in infinity as a boundary condition is with galerkin methods.

Probably using a basis like $F_k(x)=e^{-\alpha x^2}\sin(k x)$. Though you should look more into the probable basis.

Here the boundary conditions are implicit in the choise of the basis, since each function already satisfies them.

Then you should do a series expansion of your initial conditions (like in a fourier series, but with the new basis). witch will give you a system of ODE-s for the coeficients of the series expansion.

I hope this is of any help.