[Mathematica 9] Численное решение системы УрЧП

[kw_artem]

Здравствуйте! Решил реализовать на Математике решение системы уравнений:

$ \begin{cases} a_t=D_1a_{xx}+K_1a b (1-a) \frac{1+K_2 a}{1+K_3 h}-a, \\ b_t=D_2b_{xx}+a-K_4 b, \\ h_t=D_3h_{xx}+K_5 a^2-K_6h; \\ \end{cases} $

с граничными и начальными условиями:

$\displaystyle a_x(0,t)=a_x(L,t)=b_x(0,t)=b_x(L,t)=h_x(0,t)=h_x(L,t)=0,$
$a(x,0)=w(x), b(x,0)=0, h(x,0)=0.$

Написал соответствующий код:

(*=====Коэффициенты и размер области======*)
K1 = 6.85; K2 = 10; K3 = 2.36; K4 = 0.9; K5 = 20.0; K6 = 0.05;
D1 = 1; D2 = 1; D3 = 1;
L = 500; T = 240;

(*=====Функция "Шапочка" !как начальный импульс функции a(x,t)!=====*)
w[x_] := Piecewise[{{Exp[-(25/(5^2 - (x - 10)^2))], -5 < x - 10 < +5}}, 0];

s = NDSolve[
  {(*=====Система из трёх уравнений для функций a,b,h от (x,t)=====*)
   D[a[x, t], t] == 
    D1*D[a[x, t], x, x] + 
     K1*a[x, t]*
      b[x, t]*(1 - a[x, t])*((1 + K2*a[x, t])/(1 + K3*h[x, t])) - 
     a[x, t],
   D[b[x, t], t] == D2*D[b[x, t], x, x] + a[x, t] - K4*b[x, t],
   D[h[x, t], t] == D3*D[h[x, t], x, x] + K5*(a[x, t])^2 - K6*h[x, t],
   (*=====Граничные условия ∂_x=0 на концах отрезка [0,L]=====*) 
   Derivative[1, 0][a][0, t] == 0, Derivative[1, 0][a][L, t] == 0, 
   Derivative[1, 0][b][0, t] == 0, Derivative[1, 0][b][L, t] == 0, 
   Derivative[1, 0][h][0, t] == 0, Derivative[1, 0][h][L, t] == 0,
   (*b[0,t]==b[500,t],h[0,t]==h[500,t],a[0,t]==a[500,t],*)
   (*=====Начальные условия=====*)
   a[x, 0] == w[x], b[x, 0] == 0, h[x, 0] == 0},
  {a, b, h}, {x, 0, L}, {t, 0, T}]

Но Mathematica отказывается считать. Вернее она "считает" и вот, что показывает при $t=0$ для функции $a(x,t)$ в "центре" шапочки ($x=10$):

_In[18]:= a[10, 0] /.s
Out[18]:= {0.}

когда нуля там не может быть. Может я неправильно использую команду NDSolve[]? Может она заточена только под линейные УрЧП? Ведь первое и третье уравнения системы как раз нелинейны. Хотя, интересно, когда вместо написанных выше граничных условий (на производную) указал другие -- периодические условия(они приведены выше, закомментированы), код заработал. Прошу помочь разобраться. Заранее спасибо!

[Vince Diesel]

Похоже, что пик слишком локализован для такой области и значений параметров по умолчанию в Dsolve. Добавление опции StartingStepSize -> 1 помогает.

[kw_artem]

Vince Diesel, Спасибо! Действительно сработало!