Hi

I am given two coupled differential equations and try to analyze them numerically because as far as I know, there is no analytical solution known(only the steady-state on an infinite domain is known, which I use to define initial conditions).

These two equations basically describe a mass conserving reaction-diffusion system to which I may add perturbations that are geometrically induced advection terms(still mass conserving).

When I activate the perturbation I have small advection terms that shift the pattern but do not deform the pattern in a relevant way. That's also what I observe numerically.

The advection equation(without diffusion, etc) reads du/dt + v du/dx = 0

v > 0 -> pattern moves to the right.

v < 0 -> pattern moves to the left.

In my equations, the dominant advection terms are proportional to dR/dx which I called Rfirstderiv(x). So v is proportional to Rfirstderiv(x) and also its sign is determined by Rfirstderiv(x).

Also, intuitively the pattern is only sensitive to advection perturbations by Rfirstderiv(x) where the interface is located because that's where c du/dx does not vanish. So you would suspect, that the direction of the pattern shift changes when Rfirstderiv(x) changes sign at some point.

But in my numerics, the interface just moves across the extremum of R as if there wasn't any change in the sign of the advection term. This has nothing to do with the boundaries: The very same behavior also occurs when you apply periodic boundaries.

Also note that I have a two-component system where the corresponding interfaces are at the same position in x, so the argument above should be valid Here I add the corresponding notebook.

-The first cell defines the parameters, the reaction kinetics, and the stationary pattern without perturbation(coarsening neglected)

-the second cell defines the geometric perturbation and includes plots where you can see yourself, that the advection term undergoes a change in sign when crossing the extremum.

-third cell defines initial conditions

-fourth cell defines the dynamic equations and boundary conditions

-fifth cell just executes the simulation

-sixth cell shows the two components as a function of time. One of the plots also shows the geometry(radius) where the perturbation is multiplied by 1000 for visibility.

-seventh line checks mass-conservation

That kind of behavior doesn't depend on whether I use the BDF method or the IDA.

Does anyone know what's going on? Can it be correct? I suspect there is an error with the numerics but I don't know what it is.

I would be really thankful if someone could help. Thanks in advance

Luca