0
$\begingroup$

I realized that in this solution, changing one of the boundary conditions makes the answer invalid. I wonder if it is possible to modify the code in this answer to solve the following equations.

$$ \partial_t a(x,t) + \partial_x (a(x,t) u(x,t))=0 $$ $$ \partial_x (a(x,t) |u_x|^{m-1} u_x)=0$$

Initial and boundary conditions are $$ u(0,t)= exp(-t), $$

$$ a(x,0)= (L_0-x)^b , L_0=L(t=0)=5; $$ $$ a(0,t)= (L_0)^b, $$

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

$m$ is a positive real number. values of $b>0$ and $m$ are arbitrary. In particular $m$ is non-integer.

Edit: This has now been solved nicely by Alex Trounev. Unfortunately, these initial conditions lead to a not-so-interesting solution $u=u(t)$. A more interesting case would be

$𝑒(π‘₯,0)=(𝐿0βˆ’π‘₯)^{(1βˆ’π‘/π‘š)}$

and

$\partial_x 𝑒(0,𝑑)=βˆ’(1βˆ’π‘/π‘š)(𝐿0)^{βˆ’π‘/π‘š} 𝑒π‘₯𝑝(βˆ’π‘‘)$

Any ideas on how to solve this? I tried the following but I get an error.

    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;
   mu = RealAbs[u1]^(m - 1);
   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  mu  u1);
   eq2[[1]] = vU[[1]] - (L[0]^(1 - b/m))  Exp[-t];
   eq3 = {L'[t] - vU[[-1]]};
   ic = Join[vA - (L[0] - L[0]  xgrid)^b /. t -> 0, {L[0] - 5}];
   icb = 
    Join[vU - (L[0] - L[0]  xgrid)^(1 - b/m) /. t -> 0, {L[0] - 5}];
   var = Join[vA, vU, {L[t]}]; eqs = Join[eq1, eq2, eq3];
   s = NDSolve[{Table[eqs[[i]] == 0, {i, Length[eqs]}], 
      Table[ic[[i]] == 0, {i, Length[ic]}], 
      Table[icb[[i]] == 0, {i, Length[icb]}]}, var, {t, 0, tmax}, 
     Method -> {"IndexReduction" -> Automatic}]; s[[1]]];

sol11 = sol[0.5, 2];
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; xmax = x0[1];

{DensityPlot[
  If[x <= x0[t], a[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "a", PlotPoints -> 100], 
 DensityPlot[
  If[x <= x0[t], u[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "u", PlotPoints -> 100], 
 Plot[L[t] /. sol11, {t, 0, 1}, Frame -> True, 
  FrameLabel -> {"t", "L"}]}

Now, I have used new code to do this but I still get errors

    sol[b_, m_] := Module[{tmax = 1, dy = 1/10}, xgrid = Range[0, 1, dy];
   nn = Length[xgrid];
   M2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
      DifferenceOrder -> 4]@"DifferentiationMatrix";
   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;
   mu = RealAbs[u1]^(m - 1);
   eq1 = D[vA, t] - (L'[t]/L[t])  xgrid  a1 + 1/L[t]  M1 . (vA  vU);
   eq1[[1]] = D[vA[[1]], t];
   eq2 = D[M1 . (vA  mu  u1), t];
   eq2[[1]] = D[vU[[1]], t] - (1 - b/m) (L[0])^(-b/m)  Exp[-t];
   eq3 = {L'[t] - vU[[-1]]};
   (*ic=Join[vA-(L[0]-L[0] xgrid)^b/. t->0,{L[0]-5},vU-.1 xgrid/. t->
   0];var=Join[vA,vU,{L[t]}];*)
   ic = Join[vA - (L[0] - L[0]  xgrid)^b /. t -> 0, {L[0] - 5}, 
     vU - (L[0] - L[0]  xgrid)^(1 - b/m) /. t -> 0]; 
   var = Join[vA, vU, {L[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]]];
sol55 = sol[.5, 1.3];

lstA = Table[{{xgrid[[i]], t}, vA[[i]]} /. sol55, {i, 
   Length[xgrid]}, {t, 0, 1, .02}]; a = 
 Interpolation[Flatten[lstA, 1]]; lstU = 
 Table[{{xgrid[[i]], t}, vU[[i]]} /. sol55, {i, Length[xgrid]}, {t, 0,
    1, .02}]; u = Interpolation[Flatten[lstU, 1]];
x0[t_] = L[t] /. sol55; xmax = FindMaximum[x0[t], {t, .5}][[1]];
{DensityPlot[
  If[x <= x0[t], a[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "a", PlotPoints -> 100], 
 DensityPlot[
  If[x <= x0[t], u[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "u", PlotPoints -> 100], 
 Plot[L[t] /. sol55, {t, 0, 1}, Frame -> True, 
  FrameLabel -> {"t", "L"}]}
$\endgroup$
2
  • 1
    $\begingroup$ At a minimum, please provide you equations in Mathematica format. Also, are you seeking a symbolic or numerical solution? If the latter, please provide the values of all constants. $\endgroup$
    – bbgodfrey
    Commented Apr 25 at 13:56
  • $\begingroup$ I have added a comment on this. Thanks for your comment @bbgodfrey $\endgroup$ Commented Apr 25 at 14:04

1 Answer 1

1
$\begingroup$

We can solve this problem using FDM code from here with a small modification

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; 
   mu = RealAbs[u1]^(m - 1);
   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 mu u1);
   eq2[[1]] = vU[[1]] - Exp[-t]; eq3 = {L'[t] - vU[[-1]]};
   ic = Join[vA - (L[0] - L[0] xgrid)^b /. t -> 0, {L[0] - 5}]; 
   var = Join[vA, vU, {L[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 -> {"IndexReduction" -> Automatic}]; s[[1]]];

Example of usage

sol11 = sol[2, 2];

Visualization

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; xmax = x0[1];

{DensityPlot[
  If[x <= x0[t], a[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "a", PlotPoints -> 100], 
 DensityPlot[
  If[x <= x0[t], u[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "u", PlotPoints -> 100], 
 Plot[L[t] /. sol11, {t, 0, 1}, Frame -> True, 
  FrameLabel -> {"t", "L"}]}

Figure 1

Update 1. We can improve code above so that we can add initial data for u[x,t] as well. We have

sol[b_, m_] := 
  Module[{tmax = 1, dy = 1/10}, xgrid = Range[0, 1, dy]; 
   nn = Length[xgrid];
   M2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], xgrid, 
      DifferenceOrder -> 4]@"DifferentiationMatrix"; 
   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; 
   mu = RealAbs[u1]^(m - 1);
   eq1 = D[vA, t] - (L'[t]/L[t]) xgrid a1 + 1/L[t] M1 . (vA vU); 
   eq1[[-1]] = D[vA[[-1]], t];
   eq2 = D[M1 . (vA mu u1), t];
   eq2[[1]] = D[vU[[1]], t] + Exp[-t]; eq3 = {L'[t] - vU[[-1]]};
   ic = Join[vA - (L[0] - L[0] xgrid)^b /. t -> 0, {L[0] - 5}, 
     vU - .1 xgrid /. t -> 0]; var = Join[vA, vU, {L[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

sol55 = sol[.5, .5];

lstA = Table[{{xgrid[[i]], t}, vA[[i]]} /. sol55, {i, 
   Length[xgrid]}, {t, 0, 1, .02}]; a = 
 Interpolation[Flatten[lstA, 1]]; lstU = 
 Table[{{xgrid[[i]] , t}, vU[[i]]} /. sol55, {i, Length[xgrid]}, {t, 
   0, 1, .02}]; u = Interpolation[Flatten[lstU, 1]];
x0[t_] = L[t] /. sol55; xmax = FindMaximum[x0[t], {t, .5}][[1]];
{DensityPlot[
  If[x <= x0[t], a[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "a", PlotPoints -> 100], 
 DensityPlot[
  If[x <= x0[t], u[x/x0[t], t], Nothing], {x, 0, xmax}, {t, 0, 1}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> {"x", "t"}, PlotLabel -> "u", PlotPoints -> 100], 
 Plot[L[t] /. sol11, {t, 0, 1}, Frame -> True, 
  FrameLabel -> {"t", "L"}]}

Figure 2

$\endgroup$
27
  • $\begingroup$ Thanks for your help. I think there might be an issue with this method as it gives a solution in which u only depends on t and not on x! This can be seen by plotting Animate[Plot[u[x/x0[t], t], {x, 0, xmax}, PlotRange -> {{0, 16}, {0, 1}}], {t, 0, 1}] . This is probably why the code doesn't work for m<1 as in this case $1/ \partial_x u$ goes to infinity. Any idea how to solve this? $\endgroup$ Commented Apr 27 at 14:28
  • $\begingroup$ @questionerno8 No it not right. Please check picture in my answer. There are clearly seen that a depends on x and t as well. :) $\endgroup$ Commented Apr 27 at 15:56
  • $\begingroup$ Sorry, my last comment was wrong but still, the first one is correct and the solution gives u=u(t) and is independent of x. this is not very interesting as one could solve the equations analytically in this case. I think this might be resolved by using $$ \partial_x u(0,t)= exp(-t), $$. Is it possible to use this boundary condition and see if it gives u to be a function of x as well? $\endgroup$ Commented Apr 27 at 16:02
  • $\begingroup$ Equations have first order derivative of x, but why doesn't that allow us to use $\partial_x u(x_0,t) =exp(-t)$ at x_0=0? $\endgroup$ Commented Apr 27 at 16:17
  • 1
    $\begingroup$ @questionerno8 You are wrong about $F$, not about $f$, also please pay attention that $F$ nothing to do with $f$ :) $\endgroup$ Commented Apr 29 at 20:32

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.