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

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?

• > the results are quite strange... Can you clarify why you think the results are strange? At least I don't have physical intuition on this problem and the plot seems rather ok. What is the issue exactly? Sep 25, 2021 at 12:20
• @HansOlo, the two equations above are decoupled since I set A as a constant. First, the equation for $Z(t, r)$ is simply a diffusion equation, in which the initial value of a Cos function will be smoothed and finally decays to 0. So the first figure is correct. But then the integral of a Cos[r] in space should be Sin[r]. Therefore the second plot is not correct. You can found those strange bump up around $t=0$ and $r\approx 1$ in the second plot. Sep 25, 2021 at 13:12
• @YanQH It is not clear why NDSolve getting solution for Z since there is no any special algorithm to solve mixture integral and differential equations, see, for example, my post about it on mathematica.stackexchange.com/questions/217201/… Sep 26, 2021 at 7:16
• @AlexTrounev NDSolve seems like solved the system as DAEs, and satisfy those boundary condition strictly. But the initial value of $V(t, r)$ isn't what I expected. As for your post you refereed, I think there might be an error in the integral. Integrate[W[x, y, t], {y, 0, y}] doesn't make sense to me. I would write it as Integrate[W[x, yy, t], {yy, 0, y}] or Integrate[W[x, y, t], {y, 0, L}] depending on the problem. Sep 26, 2021 at 12:29
• @YanQH There is no special algorithm to solve mixture from integral and differential equations in NDSolve. But we can try to implement some method like FDM or Method of Lines for this case. See my answer with FDM algorithm. Sep 29, 2021 at 16:41

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

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

fd1 = NDSolveFiniteDifferenceDerivative[Derivative[1], grid,
DifferenceOrder -> 8]; m1 = fd1["DifferentiationMatrix"]; fd2 =
NDSolveFiniteDifferenceDerivative[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"}]
`

• Thank you very much for such detailed answer. Since I only started to use MMA 1 year ago, I still have many things to learn. Sometimes, procedures of NDSolve work like black boxes. So I am curious and think there might be a chance that it will work for this special case. Anyway, a FDM algorithm like yours could be used more universally, and also more suitable for my real problem. It's very helpful for me to see a workable example with MMA's FDM packages. I can learn to build my algorithm from it. Thank you again! Sep 30, 2021 at 16:07
• @YanQH You are welcome! Oct 1, 2021 at 6:21