0
$\begingroup$

Please, help me with the following issue. I am solving the following system of the coupled PDEs

PDEs

with the following boundary conditions

Boundary

and initial conditions.

Initial

After performing the following simplified code (for simplicity I set all parameters equal to 1):

(*specify the system*)
op = 
  {-Inactive[Grad][n[t, x], t] + Laplacian[n[t, x], x] - 
     Inactive[Grad][n[t, x]*Grad[a[t, x], x], x] + a[t, x]*b[t, x]  -  
     n[t, x] - (n[t, x])^2,
   -Inactive[Grad][a[t, x], t] + Inactive[Laplacian][a[t, x],x] + 
     (1 + Tanh[1 - b[t, x]] ) - (1 + b[t, x])*a[t, x],
   -Inactive[Grad][b[t, x], t] + 
      Inactive[Grad][n[t, x]*Inactive[Grad][b[t, x], x], x] - 
      Inactive[Grad][n[t, x], x] + n[t, x]*Inactive[Grad][a[t, x], x]}

(*specify the Dirichlet boundary conditions*)
bcs = 
  {DirichletCondition[n[t, x] == E^-t, x == 0], 
   DirichletCondition[b[t, x] == 1, x == 0]};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = {NeumannValue[-a[t, x], x == 0], NeumannValue[0, x == 1]};
BNeumn = NeumannValue[0, x == 1];
BNeumb = NeumannValue[0, x == 1];
 
(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs = 
 {n[0, x] == Piecewise[{{n0[x], x < xTild}, {0, x > xTild}}], 
  a[0, x] == 0, 
  b[0, x] == Piecewise[{{b0[x], x < xTild}, {0, x > xTild}}]};

(*solve the system*)
{nfun, afun, bfun} = 
  NDSolveValue[
    {op == {BNeumn, BNeuma, BNeumn}, bcs, incs}, 
    {n, a, b}, {x, 0, 1}, {t, 0, 10}];

I get the following error:

Set::shape: Lists {nfun, afun, bfun} and NDSolveValue[<<1>>] are not the same shape.

Please, give me a hint, why my code is incorrect. If it is nessesary, I can upload a full code with all parameters, but it gives exactly the same error. My Mathematica version is 12.1.

$\endgroup$
7
  • $\begingroup$ Looks like NDSolveValue did not return a value....(Sorry, can't check the code from a phone.) $\endgroup$
    – Michael E2
    Commented Sep 22, 2020 at 21:32
  • $\begingroup$ We need to add boundary conditions at x=1. $\endgroup$ Commented Sep 22, 2020 at 21:58
  • $\begingroup$ @MichaelE2 Good point, thanks. Are any suggestions why this happens? $\endgroup$ Commented Sep 22, 2020 at 22:17
  • $\begingroup$ @AlexTrounev Thanks, I will try. But according to the tutorial reference.wolfram.com/language/FEMDocumentation/tutorial/… if no boundary conditions are specified,NDSolveValue suggestes that they are zero Neumann (which is indeed what I want according to the equations 14 from my pictures). Or it is wrong? $\endgroup$ Commented Sep 22, 2020 at 22:19
  • 1
    $\begingroup$ Normally, I would expect there to be at least one error before the Set::shape error that would given me a hint. Incorrect syntax or failure to be able to compute the solution are common reasons. $\endgroup$
    – Michael E2
    Commented Sep 22, 2020 at 23:32

2 Answers 2

1
$\begingroup$

In this case it is better to use not FEM orientated code as follows

eq = {-D[n[t, x], t] + D[n[t, x], x, x] - 
     D[n[t, x]*D[a[t, x], x], x] + a[t, x]*b[t, x] - 
     n[t, x] - (n[t, x])^2 == 
    0, -D[a[t, x], t] + 
     D[a[t, x], x, x] + (1 + Tanh[1 - b[t, x]]) - (1 + b[t, x])*
      a[t, x] == 
    0, -D[b[t, x], t] + D[n[t, x]*D[b[t, x], x], x] - D[n[t, x], x] + 
     n[t, x]*D[a[t, x], x] == 0};

(*specify the Dirichlet boundary conditions*)
bcs = {n[t, x] == E^-t /. x -> 0, 
   b[t, x] == 1 /. 
    x -> 0, -Derivative[0, 1][a][t, x] + a[t, x] == 0 /. x -> 0, 
   Derivative[0, 1][n][t, 1] == 0, Derivative[0, 1][a][t, 1] == 0, 
   Derivative[0, 1][b][t, 1] == 0};

(*specify the Neumann boundary conditions (for function a(t,x))*)
BNeuma = NeumannValue[-a[t, x], x == 0];

(*function to use in the picewise intial conditions*)
xTild = 0.05;
n0[x_] := 1/xTild^3*(x - xTild)*(2*x^2 - xTild*x - xTild^2);
b0[x_] := 1/xTild^3 (x - xTild)*(2*x^2 - xTild*x - xTild^2);

(*initial conditions themselfs*)
incs = {n[0, x] == If[x <= xTild, n0[x], 0], a[0, x] == 0, 
   b[0, x] == If[x < xTild, b0[x], 0]};


{nfun, afun, bfun} = 
  NDSolveValue[Join[eq, bcs, incs], {n, a, b}, {x, 0, 1}, {t, 0, 10}];

There is message NDSolveValue::mxsst: Using maximum number of grid points 10000 allowed by the MaxPoints or MinStepSize options for independent variable x., but numerical solution exists and we can visualize it as

{DensityPlot[nfun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All], 
 DensityPlot[afun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All], 
 DensityPlot[bfun[t, x], {x, 0, 1}, {t, 0, 10}, 
  ColorFunction -> "Rainbow", PlotLegends -> Automatic, 
  FrameLabel -> Automatic, PlotRange -> All]}

Figure 1

$\endgroup$
4
  • 1
    $\begingroup$ Thank you very much! But how did you conclud that it is better to use not FEM? Is it clear from the system itself? $\endgroup$ Commented Sep 23, 2020 at 11:45
  • $\begingroup$ @АндрейКокорев It is too complicated to use FEM for system of nonlinear PDE. We know that for parabolic systems there is method of lines, So we can solve it without FEM. $\endgroup$ Commented Sep 23, 2020 at 12:25
  • $\begingroup$ thanks again, you helped a lot! $\endgroup$ Commented Sep 23, 2020 at 13:53
  • $\begingroup$ @АндрейКокорев You are welcome! $\endgroup$ Commented Sep 23, 2020 at 13:55
0
$\begingroup$

As was already said, NDSolve does not evaluate. The reason is that the syntax is wrong. You should specify something like:

NDsolve[{equations, boundary conditions},..]

But you specified:

NDsolve[{{expression1,Expression2}=={ boundary conditions},..]

You must give equations, not expressions.

$\endgroup$

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.