Nieliniowy system do chemotaksji

Question

Chcę rozwiązać tryb chemotaksji, podany przez następny system nieliniowy :

Pochodzi z książki Murraya: równanie (11.30) na pag. 408

$$\frac{\partial n}{\partial t} = D \frac{\partial^2 n}{\partial x^2} -\xi_0 \partial_x \Bigl( n \frac{\partial a}{\partial x} \Bigr)$$

$$\frac{\partial a}{\partial t} = hn - ka + D_a \frac{\partial^2 a}{\partial x^2}$$

gdzie $h,k,D_a,D$ to tylko parametry i $D_a>D$ a domena to $x \in [-6,6]$

Postanowiłem przyjąć warunki brzegowe bez strumienia , tj$$\partial_x(n(-6,t))=\partial_x (a(-6,t))=0$$ $$\partial_x(n(6,t))=\partial_x (a(6,t))=0$$

i jako warunki początkowe $$n(0,x)=e^{-x^2}$$ $$a(0,x)=\cos( \pi x)$$

Zauważ, że liczbowo warunki są zgodne, ponieważ wykładnia jest „płaska”. Wiem, że analitycznie to nieprawda.

Zintegrowałem się do czasu $T=0.1$ własnym solwerem MES (z liniowymi elementami skończonymi) i uzyskaj, korzystając z parametrów $$D = 2 \quad D_a = 5.5 \quad h = 0.5 \quad k = 0.5 \quad \xi_0 = 0.2$$

Chciałbym użyć Mathematica, aby sprawdzić moje wyniki i wypróbować zmiany niektórych parametrów, ale nie mogę zrozumieć, jak rozwiązać nieliniowy układ, taki jak powyższy. Czy ktoś mógłby pokazać działkę, którą powinienem uzyskać za pomocą Mathematica i, jeśli to możliwe, właściwy fragment kodu?

---

EDYTOWAĆ:

Oto, co otrzymałem, które ma kształt odpowiedzi Daniela, która wydaje się być podobna do jego

EDYTOWAĆ:

Podstawową zasadą modelu jest:

Ameby śluzowca Dictyostelium discoideum o gęstości n (x, t) wydzielają chemiczny atraktant, cykliczny AMP i zaczynają tworzyć się przestrzenne skupiska ameb. Książka mówi, że ti używa warunków brzegowych o zerowym strumieniu i to jest w porządku. Ale jakie warunki początkowe mogę wykorzystać$n(x,t)$ i $a(x,t)$ które są fizycznie istotne?

Answer 1

Jeśli używasz metody elementów skończonych, no fluxjest to domyślny warunek brzegowy, więc nie ma potrzeby określania. Alternatywą dla odpowiedzi Daniela byłoby:

(* Define parameters *)
l = 6;
tend = 0.1;
parms = {d -> 2, da -> 5.5, h -> 0.5, k -> 0.5, x0 -> 0.2};
(* Create Parametric PDE operators for n and a *)
parmnop = 
  D[n[t, x], t] - d D[n[t, x], x, x] + x0 D[n[t, x] D[a[t, x], x], x];
parmaop = D[a[t, x], t] - da D[a[t, x], x, x] + k a[t, x] - h n[t, x];
(* Setup PDE System *)
pden = (parmnop == 0) /. parms;
pdea = (parmaop == 0) /. parms;
icn = n[0, x] == Exp[-x^2];
ica = a[0, x] == Cos[π x];
(* Solve System *)
{nif, aif} = 
  NDSolveValue[{pden, pdea, icn, ica}, {n, a}, {t, 0, tend}, {x, -l, 
    l}, Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"FiniteElement", 
       "MeshOptions" -> MaxCellMeasure -> 0.1}}];

(* Display results *)
Manipulate[
 Plot[{nif[t, x], aif[t, x]}, {x, -l, l}, PlotRange -> All], {t, 0, 
  tend}, ControlPlacement -> Top]

Answer 2

Oto mój kod. Niestety, przy t == 0,1 nie powiela wyniku. Mam nadzieję, że się nie pomyliłem.

eq = {D[n[x, t], t] == 
     d  D[n[x, t], {x, 2}] - c0 D[n[x, t] D[a[x, t], x], x],
    D[a[x, t], t] == h  n[x, t ] - k a[x, t] + da  D[a[x, t], {x, 2}],
    (D[n[x, t], x] /. x -> -6) == 0, (D[a[x, t], x] /. x -> -6) == 
     0, (D[n[x, t], x] /. x ->   6) == 
     0, (D[a[x, t], x] /. x ->   6) == 0,
    n[x, 0] == Exp[-x^2], a[x, 0] == Cos[Pi x]} /. {d -> 2, da -> 5.5,
     h -> 0.5, k -> 0.5, c0 -> 0.2};
sol[x_] = {n[x, 0.1], a[x, 0.1]} /. 
  NDSolve[eq, {n, a}, {t, 0, 0.1}, {x, -6, 6}][[1]]
Plot[sol[x], {x, -6, 6}, PlotRange -> All]