ระบบ NonLinear สำหรับ chemotaxis

Question

ฉันต้องการแก้โหมด chemotaxis ซึ่งกำหนดโดยระบบที่ไม่ใช่เชิงเส้นถัดไป:

นำมาจากหนังสือของ Murray: สมการ (11.30) ที่หน้า 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}$$

ที่ไหน $h,k,D_a,D$ เป็นเพียงพารามิเตอร์และ $D_a>D$ และโดเมนคือ $x \in [-6,6]$

ฉันตัดสินใจที่จะไม่ใช้เงื่อนไขขอบเขตของฟลักซ์กล่าวคือ$$\partial_x(n(-6,t))=\partial_x (a(-6,t))=0$$ $$\partial_x(n(6,t))=\partial_x (a(6,t))=0$$

และเป็นเงื่อนไขเริ่มต้น $$n(0,x)=e^{-x^2}$$ $$a(0,x)=\cos( \pi x)$$

โปรดทราบว่าเงื่อนไขเป็นตัวเลขที่เข้ากันได้เนื่องจากเลขชี้กำลังเป็น "แบน" ฉันรู้ว่าในทางวิเคราะห์แล้วมันไม่เป็นความจริง

ฉันบูรณาการได้ตลอดเวลา $T=0.1$ ด้วยตัวแก้ FEM ของฉันเอง (ที่มีองค์ประกอบ จำกัด เชิงเส้น) และรับสิ่งต่อไปนี้โดยใช้พารามิเตอร์ $$D = 2 \quad D_a = 5.5 \quad h = 0.5 \quad k = 0.5 \quad \xi_0 = 0.2$$

ฉันต้องการใช้ Mathematica เพื่อตรวจสอบผลลัพธ์ของฉันและลองสิ่งที่ออกมาโดยการเปลี่ยนพารามิเตอร์บางอย่าง แต่ฉันไม่เข้าใจวิธีแก้ปัญหาระบบที่ไม่ใช่เชิงเส้นเช่นเดียวกับด้านบน มีใครสามารถแสดงพล็อตที่ฉันควรได้รับจาก Mathematica และถ้าเป็นไปได้ข้อมูลโค้ดที่ถูกต้อง

---

แก้ไข:

นี่คือสิ่งที่ฉันได้รับซึ่งมีรูปร่างของคำตอบของแดเนียลซึ่งดูเหมือนจะคล้ายกับของเขา

แก้ไข:

หลักการทางกายภาพที่อยู่เบื้องหลังแบบจำลองคือ:

อะมีบาของราเมือก Dictyostelium discoideum ที่มีความหนาแน่น n (x, t) หลั่งสารดึงดูดทางเคมีไซคลิกแอมป์และการรวมตัวเชิงพื้นที่ของอะมีบาเริ่มก่อตัว หนังสือเล่มนี้บอกว่าใช้เงื่อนไขขอบเขตศูนย์ฟลักซ์และไม่เป็นไร แต่ฉันสามารถใช้เงื่อนไขเริ่มต้นอะไรได้$n(x,t)$ และ $a(x,t)$ ที่เกี่ยวข้องกับร่างกาย?

Answer 1

หากคุณใช้วิธีไฟไนต์เอลิเมนต์no fluxเป็นเงื่อนไขขอบเขตเริ่มต้นดังนั้นจึงไม่จำเป็นต้องระบุ ทางเลือกอื่นสำหรับคำตอบของ Daniel คือ:

(* 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

นี่คือรหัสของฉัน น่าเสียดายที่ t == 0.1 ผลลัพธ์ของคุณไม่ซ้ำกัน ฉันหวังว่าฉันจะไม่ทำผิดพลาด

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]