0
$\begingroup$

I am trying to solve the following equations in the interval $0<t<10$ and $0<y<1$.

$\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?

enter link description here

$\endgroup$
3

1 Answer 1

3
$\begingroup$

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 = NDSolve`FiniteDifferenceDerivative[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"}]

Figure 1

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}]}

Figure 2

$\endgroup$
8
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jul 10 at 10:59
  • $\begingroup$ @questionerno8 Yes you are right. I have revised code and added new pictures. $\endgroup$ Commented Jul 10 at 11:18
  • $\begingroup$ 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. $\endgroup$ Commented Jul 10 at 11:54
  • $\begingroup$ 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? $\endgroup$ Commented Jul 10 at 12:19
  • $\begingroup$ I agree that it is OK to use two boundary conditions on one boundary. $\endgroup$ Commented Jul 10 at 15:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.