error with NDSolve of PDE with Neumann conditions

Question

I'm trying to solve the nonlinear "reaction-diffusion" system with Neumann conditions:

Subscript[r,a]=3;Subscript[r,h]=3/2;\[Delta]=2;Subscript[\[Alpha],h]=2;Subscript[\[Alpha],a]=2;
    eqns40f={Subscript[\[Alpha], a] D[A[w,t],{w,2}]+Subscript[\[Alpha],a]^2/Subscript[\[Alpha], h] A[w,t]^2/H[w,t]-Subscript[r, a]A[w,t]+Subscript[r, a]==D[A[w,t],t],\[Delta] Subscript[\[Alpha], h] D[H[w,t],{w,2}]+\!\(
    \*SubsuperscriptBox[\(\[Alpha]\), \(a\), \(2\)]\ 
    \*SuperscriptBox[\(A[w, t]\), \(2\)]\)-Subscript[r, h] H[w,t] +Subscript[r, h]==D[H[w,t],t],Derivative[1,0][A][0,t]==0,Derivative[1,0][A][\[Pi],t]==0,Derivative[1,0][H][0,t]==0,Derivative[1,0][H][\[Pi],t]==0,A[w,0]==Cos[2w],H[w,0]==Cos[3 w]}/.{Subscript[\[Alpha], h]-> 2, Subscript[\[Alpha], a]->3}
    sol40f=NDSolve[eqns40f,{A,H},{w,0,\[Pi]},{t,0,100}];

$r_a=3;\alpha _h=2;r_h=\frac{3}{2}; \delta = 2,\alpha _a = 3;$ $\text{eqns40f}=\left\{\frac{\alpha _a^2 A(w,t)^2}{\alpha _h H(w,t)}-r_a A(w,t)+\alpha _a \frac{\partial ^2A(w,t)}{\partial w^2}+r_a=\frac{\partial A(w,t)}{\partial t},\alpha _a^2 A(w,t)^2-r_h H(w,t)+\delta \alpha _h \frac{\partial ^2H(w,t)}{\partial w^2}+r_h=\frac{\partial H(w,t)}{\partial t},A^{(1,0)}(0,t)=0,A^{(1,0)}(\pi ,t)=0,H^{(1,0)}(0,t)=0,H^{(1,0)}(\pi ,t)=0,A(w,0)=\cos (2 w),H(w,0)=\cos (3 w)\right\};$ $\text{sol40f}=\text{NDSolve}[\text{eqns40f},\{A,H\},\{w,0,\pi\},\{t,0,100\}];$

But the next errors emerge:

NDSolve::ibcinc: Warning: boundary and initial conditions are inconsistent.

NDSolve::ndsz: At t == 1.047627475220518`*^-17, step size is effectively zero; >singularity or stiff system suspected.

NDSolve::eerr: Warning: scaled local spatial error estimate of 1163.0196390315273at t = 1.047627475220518*^-17 in the direction of independent variable w is much greater than the prescribed error tolerance. Grid spacing with 29 points may be too large to achieve the desired accuracy or precision. A singularity may have formed or a smaller grid spacing can be specified using the MaxStepSize or MinPoints method options.

Boundary and initial conditions are actually consistent.

Answer 1

The main problem here, is that H is in the denominator in one of the terms of the pde, and H is zero at certain values of w in the initial condition. The result is that NDSolve quits at t = 10^-17 which is essentialy t = 0 for machine precision. I have not seen NDSolve succeed in such circumstances. We can do something, however, if there is any interest in later time behavior of A and H. Rewriting the problem without subscripts and using using clearer symbols for this forum:

ra = 3;
rh = 3/2;
delta = 2;
alphah = 2;
alphaa = 3;

pde = {
   alphaa D[A[w, t], w, w] + (alphaa^2 A[w, t]^2)/(alphah H[w, t]) - 
     ra A[w, t] + ra == D[A[w, t], t],
   delta alphah D[H[w, t], w, w] + alphaa^2 A[w, t]^2 - rh H[w, t] + 
     rh == D[H[w, t], t]
   };

bc = {
   (D[A[w, t], w] /. w -> 0) == 0,
   (D[A[w, t], w] /. w -> Pi) == 0,
   (D[H[w, t], w] /. w -> 0) == 0,
   (D[H[w, t], w] /. w -> Pi) == 0
   };

Now change the ic so that H starts positive everywhere. The shape of H is unchanged.

ic = {
  A[w, 0] == Cos[2 w],
  H[w, 0] == 1.0000001 + Cos[3 w]
  }

sol = NDSolve[{pde, bc, ic}, {A, H}, {w, 0, Pi}, {t, 0, 100}, 
  PrecisionGoal -> Infinity, MaxStepFraction -> 1/200, 
  MaxSteps -> {50000, Automatic}];

A[w_, t_] = A[w, t] /. sol[[1]];
H[w_, t_] = H[w, t] /. sol[[1]];

Plot the solution.

gifs = Table[
   Plot[{A[w, t], H[w, t]}, {w, 0, Pi}, PlotRange -> {-3, 13}, 
    PlotLegends -> {"A", "H"}], {t, 0, 4.1, .05}];
ListAnimate[%]

What is interesting, if we now change the ic on H to be negative everywhere, i.e.

ic = {
  A[w, 0] == Cos[2 w],
  H[w, 0] == -1.0000001 + Cos[3 w]
  }

we get:

Crossing zero in this second case seems to cause no trouble, and the two cases appear to end up with the same values after sufficient time. Fwiw, starting H at a higher center value such as 20, causes it to fall below the steady state and then rise to the same value as the two cases plotted here.

If we change the initial shape of the A and H curves, their behavior is substantially different from the cases plotted here.