# Numerical solution with both Neuman and Drichlet BC at one point

I am trying to solve the following equations in the interval $$0 and $$0.

$$\partial_t a(y,t) + \partial_y[ a(y,t) u(y,t)]=0,$$

$$\partial_y[ a(y,t) \: (m+ 3 \frac{\partial_y u(y,t)}{L(t)}+\frac{3 }{2}\frac{\partial_y \partial_t u(y,t)}{L(t)})]=0,$$

$$\partial_t L(t)= u(1,t),$$

where $$m$$ and $$b$$ are arbitrary constants.

The boundary conditions and the initial conditions are

$$\partial_y u(0,t)=0,$$

$$u(0,t)=1,$$

$$L(t=0)=5,$$

$$a(0,t)=5^b$$

$$a(y,t=0)=(5-y)^b$$

$$u(y,0)=1+y/20,$$

In the following link, a very similar question has been answered nicely. The difference here is in the boundary conditions. More specifically, I do not know how to add the derivative boundary condition $$\partial_y u(0,t)=0$$. Could someone help me with this?

• Link you provided is wrong. Do you mean this one mathematica.stackexchange.com/questions/304012/… ? Commented Jul 10 at 2:35
• I corrected the link, thanks for pointing out! @AlexTrounev Commented Jul 10 at 2:43

Using my code from here we have

ClearAll["*"]

sol[b_, m_] := Module[{tmax = 1, dy = 1/10}, xgrid = Range[0, 1, dy];
nn = Length[xgrid];
M1 = NDSolveFiniteDifferenceDerivative[Derivative[1], xgrid,
DifferenceOrder -> 4]@"DifferentiationMatrix";
vA = Table[va[i][t], {i, nn}]; a1 = M1 . vA;
vU = Table[vu[i][t], {i, nn}]; u1 = M1 . vU;
eq1 = D[vA, t] - (L'[t]/L[t])  xgrid  a1 + 1/L[t]  M1 . (vA  vU);
eq1[[1]] = D[vA[[1]], t];
eq2 = M1 . (vA  (m + 3  u1/L[t] + 3  D[u1, t]/(2  L[t])));
eq3 = {L'[t] - vU[[-1]]};
eq2[[2]] = u1[[1]];
eq2[[1]] = D[vU[[1]], t];
ic = Join[vA - (L[t] - L[t]  xgrid)^b, {L[0] - 5},
vU - (1 + xgrid/20)] /. t -> 0;
var = Join[vA, vU, {L[t]}]; var1 = D[var, t];
eqs = Join[eq1, eq2, eq3];
s = NDSolve[{Table[eqs[[i]] == 0, {i, Length[eqs]}],
Table[ic[[i]] == 0, {i, Length[ic]}]}, var, {t, 0, tmax},
Method -> {"EquationSimplification" -> "Residual"}]; s[[1]]];


Example of usage

sol11 = sol[2, .75];


Visualization

Plot[L[t] /. sol11, {t, 0, 1}, AxesLabel -> {"t", "L"}]


Animation

lstA = Table[{{xgrid[[i]], t}, vA[[i]]} /. sol11, {i,
Length[xgrid]}, {t, 0, 1, .02}]; a =
Interpolation[Flatten[lstA, 1]]; lstU =
Table[{{xgrid[[i]], t}, vU[[i]]} /. sol11, {i, Length[xgrid]}, {t, 0,
1, .02}]; u = Interpolation[Flatten[lstU, 1]];
x0[t_] = L[t] /. sol11;
{Animate[
Plot[Sqrt[a[x/x0[t], t]], {x, 0, x0[t]},
PlotRange -> {{0, 10}, {0, 2}}], {t, 0, 1}],
Animate[Plot[u[x/x0[t], t], {x, 0, x0[t]},
PlotRange -> {{0, 10}, {-2, 2}}], {t, 0, 1}]}


• I think there might be a mistake in your code. the boundary condition is $\partial_y u=0$. You have set $\partial_t \partial_y u=0$. These are different things as $\partial_y u$ is not zero at the initial condition so with your code it stays non-zero. Commented Jul 10 at 10:59
• @questionerno8 Yes you are right. I have revised code and added new pictures. Commented Jul 10 at 11:18
• Also, since you are putting this new boundary condition on the second element of the second equation, it means that the second element of equation 2 is not satisfied anymore. Commented Jul 10 at 11:54
• It is ok if we use two elements on one border and nothing on another. On the other hand, I don't understand why you try to use 2 boundary conditions one border and nothing on another? Commented Jul 10 at 12:19
• I agree that it is OK to use two boundary conditions on one boundary. Commented Jul 10 at 15:24