multiplying equation by "x" changes the solution

[wd15]

The solution to the following changes when the equation is multiplied through by "x"

from fipy import Grid1D, CellVariable
from scipy import exp
from fipy import TransientTerm, ConvectionTerm, ImplicitSourceTerm, Viewer

xmax = 5.
dx = 0.05
nx = xmax/dx
mesh = Grid1D(dx=dx, Lx=xmax)

x = mesh.getCellCenters()[0]
X = mesh.getFaceCenters()[0]

numerical = CellVariable(name="numerical", mesh=mesh)

numerical.constrain(1, mesh.facesLeft)
numerical.constrain(0, mesh.facesRight)
numerical.faceGrad.constrain([0], mesh.facesLeft)

analytical = CellVariable(name="analytical value", mesh=mesh)
analytical.setValue(exp(-x**2/2.))

eqn = TransientTerm(coeff=1 / x) == \
  - ConvectionTerm(coeff=[[1]]) \
  - ImplicitSourceTerm(coeff=x)

numerical.setValue(exp(-x**2/2.))

vi = Viewer(vars=(numerical,analytical),datamin=0,datamax=2)

for i in range(10000):
eqn.solve(var=numerical,dt=1.)
vi.plot()
raw_input('stopped')

Imported from trac ticket #466, created by wd15 on 09-20-2012 at 12:27, last modified: 09-21-2012 at 17:43

[wd15]

The reason why this doesn't work is due to the pathological nature of
the equation rather than any numerical issues. Obviously, the boundary
condition on the LHS of the domain can have no influence as the
equation is simply $\phi_t=0$ at $x=0$. Assuming we start with
$\phi=1$ as the initial condition, any drift away from $1$ can't be
prevented. There are no boundary conditions on the LHS of the domain,
only the initial condition and, hence, the only thing preventing the
drift away from $1$ on the LHS of the domain is the small value on the
RHS of the equation. Unless the cell is exactly at $x=0$ then drift
will occur. The drift reduces like the size of the grid. Technically
an infinitely fine grid would solve this problem.

Trac comment by wd15 on 09-21-2012 at 17:43