Integral in NDSolve, boundary conditions and initial conditions can't both satisfied

Question

I have serval reaction-diffusion equations to be solved. Which are equivalent to the following problem:

First is a simple diffusion equation of $Z(t,r)$ like below: $$ \dfrac{\partial Z}{\partial t}=A(t,r)\dfrac{\partial^2 Z}{\partial r^2} $$

Then $Z(t,r)$ is actually defined from: $$ Z(t,r) \equiv \dfrac{\partial V(t,r)}{\partial r}$$ In my problem, I must have the value of $V(t,r)$ at every $(t,r)$ to determine the value of, for example, $A(t,r)=F(V(t,r))$. Bondary condtions are $$ V(t,0)=0,V(t,1)=0,\\ Z_r(t,0)=0,Z_r(t,1)=0 $$

For testing, I just set $A$ a constant, and tried something like below:

A = 0.01; L = 1; tmax = 10;

eqns = {V[t, r] == Integrate[Z[t, r], r],
        D[Z[t, r], t] - A D[Z[t, r], {r, 2}] == NeumannValue[0, r == 0 || r == 1]};
INIs = {Z[0, r] == 2 Pi Cos[2 Pi r]};
BCs = {DirichletCondition[V[t, r] == 0, r == 0 || r == 1 ]};
sys = {eqns, BCs, INIs};

sol = NDSolve[sys, {V, Z}, {r, 0, L}, {t, 0, tmax}]

The codes run successfully without warning, but the results are quite strange: In the result $Z(t, r)$ is correctly calculated, while $V(t, r)$ is not correct. Does anyone know why it behaviors like this? Or any other way to do deal with this problem? Thanks a lot.

Update

I found that this system is solved by NDSovle as differential algebraic equations. I tried index reduction of DAE as described in this. Here are the codes:

A = 0.01; L = 1; tmax = 10;

eqns = {V[t, r] == Integrate[Z[t, r], r],       
   D[Z[t, r], t] - A D[Z[t, r], {r, 2}] == 
    NeumannValue[0, r == 0 || r == L]};
INIs = {Z[0, r] == 2 Pi Cos[2 Pi r],
    V[0, r] == Sin[2 Pi r]};
BCs = {DirichletCondition[V[t, r] == 0, r == 0 || r == L]
   };
sys = {eqns, BCs, INIs};

sol = NDSolve[sys, {V, Z}, {r, 0, L}, {t, 0, tmax}, 
  Method -> {"IndexReduction" -> {"StructuralMatrix"}}]

And the results are different.

• In the previous codes, all BCs are satisfied, but initial value of $V(t,r)$ is not I wanted.
• In the updated codes, with index reduction. I can add the initial value of $V(t,r)$, which is satisfied by the result. But one of the boundary conditions of $Z(t,r)$ is violated.

Why can't it satisfy both ICs and BCs as I gave?

Answer 1

We can recommend FDM algorithm of 8 order with using DAE solver as follows

Clear["Global`*"]
n = 137; grid = Range[0, n]/n;

fd1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], grid, 
  DifferenceOrder -> 8]; m1 = fd1["DifferentiationMatrix"]; fd2 = 
 NDSolve`FiniteDifferenceDerivative[Derivative[2], grid, 
  DifferenceOrder -> 8]; m2 = fd2["DifferentiationMatrix"]; varz = 
 Table[z[i][t], {i, Length[grid]}]; varz1 = 
 Table[z[i]'[t], {i, Length[grid]}]; zxx = m2 . varz; varZ = 
 Table[z[i], {i, Length[grid]}]; varv = 
 Table[v[i][t], {i, Length[grid]}]; vx = m1 . varv; varV = 
 Table[v[i], {i, 2, n}];
A = .01; v[1][t_] := 0; 
v[n + 1][t_] := 0; eq1 = 
 Table[z[i]'[t] - A zxx[[i]] == 0, {i, 2, n, 1}]; eq2 = 
 Table[z[i][t] - vx[[i]] == 0, {i, 2, n, 1}]; ic = 
 Table[z[i][0] == 2 Pi Cos[2 Pi grid[[i]]], {i, 2, n}]; ic1 = 
 Table[v[i][0] == Sin[2 Pi grid[[i]]], {i, 2, n}]; bc = {z[1][t] == 
   z[2][t], z[n + 1][t] == z[n][t]};


eqn = Join[eq1, eq2, ic, ic1, bc]; var = Join[varV, varZ];

sol = NDSolve[eqn, var, {t, 0, 10}, 
   Method -> {"EquationSimplification" -> "Residual"}];

Visualization

lst = Flatten[
   Table[{t, grid[[i]], z[i][t] /. sol[[1]]}, {t, 0, 10, .1}, {i, 
     Length[grid]}], 1];

ListPlot3D[lst, Mesh -> None, ColorFunction -> Hue, PlotRange -> All, 
 AxesLabel -> {"t", "x", "Z"}]

lst1 = Flatten[
   Table[{t, grid[[i]], v[i][t] /. sol[[1]]}, {t, 0, 10, .1}, {i, 
     Length[grid]}], 1];

ListPlot3D[lst1, Mesh -> None, ColorFunction -> Hue, PlotRange -> All,
  AxesLabel -> {"t", "x", "V"}]