Równania reakcji-dyfuzji

Dotychczas analizowaliśmy modele dynamiki populacyjnej oraz kinetykę reakcji chemicznych bez uwzględnienia zjawisk rozprzestrzeniania się populacji lub substancji chemicznych na powierzchni czy w przestrzeni. Interesowało nas tylko zachowanie się układu w czasie, a nie w przestrzeni. Oznacza to, że modelowanie oparte było na równaniach różniczkowych zwyczajnych (występują tylko pochodne ze względu na czas) . Uwzględnienie zjawisk związanych z przestrzennymi zmianami modeluje się z pomocą równań różniczkowych cząstkowych (występują pochodne ze względu na zmienne przestrzenne). Poniższy przykład przybliży nam sposób takiego modelowania.

Przykład: zjawisko dyfuzji

Wrzućmy kroplę farby akwarelowej do pojemnika z wodą. Zaobserwujemy, że z upływem czasu cząsteczki farby rozpływają się w coraz to większym obszarze i po odpowiednio długim czasie wypełniaja całkowicie pojemnik z wodą. To, co obserwujemy jest zjawiskiem dyfuzji lub, jak niektórzy nazywają, ruchem cząstek Browna. Jest to zjawisko natury molekularnej, spowodowane nieregularnymi (chaotycznymi) zderzeniami cząsteczek wody i cząsteczek farby. Zauważamy też, że podwyższenie temperatury wody powoduje szybsze mieszanie się obu substancji. W języku fizyki statystycznej oznacza to, że współczynnik dyfuzji ulega zwiększeniu, rosnie natężenie fluktuacji termicznych i następuje przyspieszenie procesu dyfuzji.W języki matematyki, proces ten jest opisany za pomocą równania:

$$\frac{\partial c(\vec r, t)}{\partial t} = D \Delta c(\vec r, t) = D \left[\frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} + \frac{\partial ^2}{\partial z^2} \right] c(\vec r, t)$$

gdzie $c(\vec r, t)$ jest koncentracją substancji w chwili czasu $t$ w punkcie przestrzeni określonej przez wektor wodzący $\vec r = (x, y, z)$. Zwracamy uwagę, że proces dyfuzji jest opisany operatorem Laplace’a (dwukrotnym różniczkowaniem).

Dla początkowej koncentracji $c(x, y, z, 0)$ skupionej punktowo w początku układu odniesienia $XYZ$, rozwiązanie równania dyfuzji ma postać

$$c(x, y, z, t) = [4\pi Dt]^{-3/2} \; exp\left[-\frac{x^2+y^2+z^2}{4Dt}\right]$$

Gdyby proces przebiegał tylko w jednym kierunku przestrzennym, np. wzdłuż osi $OX$ to odpowiednie równanie ma postać

$$\frac{\partial c(x, t)}{\partial t} = D \frac{\partial^2 }{\partial x^2} c(x, t) = D \frac{\partial^2 c(x, t) }{\partial x^2}$$

Jego rozwiązaniem jest funkcja

$$c(x, t) = [4\pi Dt]^{-1/2} \; exp\left[-\frac{x^2}{4Dt}\right]$$

Otrzymujemy więc następującą receptę na uwzględnienie dyfuzyjnego (losowego) rozprzestrzeniania się populacji lub substancji chemicznej:

Do równań opisujących dynamikę populacji dodaj jeszcze wyraz z pochodnymi drugiego rzędu względem zmiennych przestrzennych. Jeżeli zmiany mogą zachodzić tylko w jednym kierunku przestrzennym, to dodajemy wyraz analogiczny jak po prawej stronie powyższego równania. Gdyby zmiany mogły zachodzić na płaszczyżnie, to dodajemy wyraz typu:

$$+ D \left[ \frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} \right] c(x, y, t)$$

Rozwiązaniem 2-wymiarowego równania dyfuzji jest funkcja

$$c(x, y, t) = [4\pi Dt]^{-2/2} \; \exp\left[-\frac{x^2+y^2}{4Dt}\right] =\frac{1}{4\pi Dt} \; \exp\left[-\frac{x^2+y^2}{4Dt}\right]$$

xxxxxxxxxx

 

1

sage: var('x,y,t')

2

sage: c(x,y,t)=1/(4*pi*t)*exp(-(x^2+y^2)/(4*t) )

3

sage: plts = [ plot3d(c(x,y,t), (x,-4,4), (y,0,8),rgbcolor = hue(t/2.0)) for t in [0.5,1.0,1.5,2.0]]

4

sage: @interact

5

sage: def _(viewer=['tachyon','jmol','canvas3d']):

6

...       show(sum(plts),viewer=viewer)

Evaluate

Language: Sage Gap GP HTML Macaulay2 Maxima Octave Python R Singular

---

Messages

sage: def plt_t(t):
...       cmax = c(0,0,t)
...       return contour_plot( c(x,y,t) ,(x,-5,5),(y,-5,5),contours=srange(0.01,cmax,cmax/20.0),cmap='spectral',figsize=(3,3))

sage: times = [.1,0.4,0.9]
sage: html.table([["t=%0.2f"%t for t in times],[plt_t(t) for t in times]])

Wszystkie dotychczas rozważane modele (Malthusa, Verhuslta, autokatalizy, Lotki-Volterry, Maya, Biełousowa-Żabotyńskiego) mogą być uogólnione na przypadek przestrzennych zmian. Dla przykładu, jeżeli równanie dynamiki populacyjnej ma postać

$$\frac{dN(t)}{dt} = F(N(t))$$

to odpowiednie równanie uwzględniające proces dyfuzji ma postać

$$\frac{\partial N(\vec r, t)}{\partial t} = F(N(\vec r, t) ) + D \Delta N(\vec r, t)$$

gdzie w zależności od zagadnienia, laplasjan jest w jednym, dwóch lub trzech wymiarach. Tego typu równanie nazywa się równaniem reakcji-dyfuzji.

Rozważamy poniżej dwa modele: model Malthusa i model Verhulsta.

Model Malthusa z migracją: Równanie Skellama

Standardowy model Malthusa to najprostszy model wzrostu lub zaniku populacji opisany przez równanie różniczkowe

$$\frac{dN(t)}{dt} = k N(t)$$

Jeżeli populacja może przemieszczać się losowo na płaszczyźnie $XY$ to ogólniona wersja tego równania ma postać

$$\frac{dN(x, y, t)}{dt} = k N(x, y, t) + D \left[ \frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} \right] N(x, y, t)$$

gdzie teraz $N(x, y, t)$ ma interpretację koncentracji populacji (substancji chemicznej) w chwili czasu $t$ w punkcie o współrzędnych $(x, y)$, czyli liczby osobników na jednostke powierzchni. Jest to jeden z prostszych przykładów układu, który nazywamy równaniem reakcji-dyfuzji. W literaturze, powyższe równanie nazywa się równaniem Skellama, który w 1951 roku zaproponował ten model jako model inwazji gatunku. Zastosował go do opisu inwazji piżmaka w Europie.

http://pl.wikipedia.org/w/index.php?title=Plik:Bisamratte-drawing.jpg&filetimestamp=20041114193425

Do Europy piżmaki sprowadzone zostały na początku XX w. z Ameryki Północnej przez księcia Colloredo-Mansfelda jako zwierzę futerkowe. Cała dzisiejsza europejska populacja dziko żyjących piżmaków pochodzi od kilku osobników, które w 1905 roku uciekły z farmy położonej w Czechach, 40 km na południowy zachód od Pragi. Na ziemiach polskich pojawił się w latach 30. XX w.

W modelu rozpatrujemy tylko takie przypadki że $N(x, y, t) \ge 0$. Nie jest to jeszcze kompletne sformułowanie problemu. Należy sprecyzować:

1. warunek początkowy $N(x, y, 0)= N_0(x, y)$

2. warunki brzegowe

Właściwe warunki brzegowe są często trudne do zdefiniowania. Na pewno na brzegach obszaru koncentracja $N(\pm \infty, t)$ powinna przyjmować ograniczone wartości. Ale to niekoniecznie musi wystarczać, aby rozwiązania posiadały poprawne własności i były pozbawione artefaktów.

Powyższe równanie reakcji-dyfuzji jest tak trudne (łatwe) do analizy jak równanie dyfuzji. Aby to zobaczyć, dokonamy podstawienia:

$$N(x, y, t)=e^{kt} c(x, y, t)$$

Wstawiając je do równania reakcji-dyfuzji, widzimy że nowa funkcja $c(x, y, t)$ spełnia równanie

$$\frac{d c(x, y, t)}{dt} = D \left[ \frac{\partial^2 }{\partial x^2} + \frac{\partial^2 }{\partial y^2} \right] c(x, y, t)$$

czyli równanie dyfuzji w dwu-wymiarowej przestrzeni. Rozwiązanie tego równania jest przedstawione powyżej. Więc funkcja $N(x, y, t)$ ma postać

$$N(x, y, t)=e^{kt} c(x, y, t) = \frac{N_0}{4\pi Dt} \; exp\left[kt -\frac{x^2+y^2}{4Dt}\right]$$

gdzie stała $N_0=N(0)$ jest liczbą osobników w populacji w chwili początkowej $t=0$.

Ponieważ koncentracja $N(x, y, t)$ populacji ma interpretację liczby osobników na jednostkę powierzchni, to

$$\int_{-\infty}^{\;\infty} \;\int_{-\infty}^{\;\infty} N(x, y, t) \,dx \,dy = N(t) = N(0) e^{kt}$$

Zauważmy, że ewolucja koncentracji populacji zależy tylko od położenia w taki oto sposób

$$N(x, y, t)= N(r ,t) = \frac{N_0}{4\pi Dt} \; exp\left[kt -\frac{r^2}{4Dt}\right]$$

gdzie

$$r^2=x^2+y^2$$

jest miarą odległości od początku układu odniesienia.

Analizę własności $N(r, t)$ można łatwo przeprowadzić korzystając z możliwości Sage. Są one przedstawione graficznie poniżej.

xxxxxxxxxx

 

1

sage: g(x,t) = (1/(4*pi*t))*exp(t-x^2/(4*t))

2

sage: @interact

3

sage: def _(t=slider(0.1,5,0.01)):

4

...       pr0 = plot( g(x,0.1),(x,0,10),color='black' )

5

...       prt = plot( g(x,t),(x,0,10),fill=True )

6

...       (pr0+prt).show(figsize=5)

Evaluate

Language: Sage Gap GP HTML Macaulay2 Maxima Octave Python R Singular

---

Messages

Obserwujemy następującą ewolucję populacji:

1. W pewnej umownej chwili początkowej rozkład populacji ma określoną formę (powyżej: dla t=0.1).

2. Wraz z upływem czasu populacja rozprzestrzenia się i jej koncentracja wokół $r=0$ (w otoczeniu wyjściowego największego skupiska) zaczyna maleć: osobniki rozchodzą się, tempo urodzin jest na razie niewielkie.

3. Po pewnym czasie koncentracja wokół $r=0$ przyjmuje minimalną wielkość, a następnie wraz z upływem czasu zaczyna wzrastać (narasta tempo urodzin). Jednocześnie populacja rozprzestrzenia się zajmując coraz to większy obszar.

4. Następnie obserwujemy, że istnieje taka chwila czasu, gdy koncentracja wokół położenia $r=0$ przekracza wartość początkową i dalej wzrasta w coraz to większym tempie, jednocześnie rozprzestrzeniając się na coraz to większym obszarze. Obserwujemy inwazję populacji. Dalej jest już jak w kiepskim horrorze. Tym nie mniej, scenariusz ten wydaje się dość rozsądny i w pewnych warunkach całkiem dobrze opisuje w ograniczonym przedziale czasowym inwazję populacji.

Pełniejsza analiza tego modelu jest przedstawiona w książce:

Nanako Shigesada and Kohkichi Kawasaki, Biological Invasion: Theory and Practice (Oxford University Press, 2001)

Model Verhulsta z migracją: Równanie Fishera-Kołmogorowa

Model Verhulsta opisuje ewolucję populacji w przypadku ograniczonych zasobów pożywienia w środowisku, w którym żyje dana populacja:

$$\frac{dn(t)}{dt} = r n(t) \; [1- n(t)]$$

gdzie przeskalowana koncentracja $n(t) = N(t)/K$ i parametr $K$ nazywa sie pojemnością środowiska.

Model ten opisuje też kinetykę dwóch reakcji autokatalitycznych:

$$A + X {\Longleftrightarrow} 2X , \quad \quad \mbox{lub reakcji} \quad A+X {\Longleftrightarrow} 2X, \quad B+X \rightarrow C$$

Jeżeli uwzględnimy proces rozprzestrzeniania się dyfuzyjnego, model ten przyjmuje postać

$$\frac{\partial n(\vec r, t)}{\partial t} = r n(\vec r, t) \;[1- n(\vec r, t)] + D \Delta n(\vec r, t)$$

W przypadku dynamiki populacyjnej, powinniśmy rozważać migrację na płaszczyżnie (laplasjan w dwóch wymiarach), a w przypadku substancji chemicznych - dyfuzję przestrzenną (pełny trójwymiarowy laplasjan). Oba realistyczne przypadki są wyjątkowo skomplikowane. Dlatego upraszczamy model i redukujemy zagadnienia do przypadku jednowymiarowej dyfuzji i poniżej rozważamy następujące równanie

$$\frac{\partial n(x, t)}{\partial t} = r n(x, t) \;[1- n(x, t)] + D \frac{\partial^2 n(x, t)}{\partial x^2}$$

Równanie to nazywa się równaniem Fishera-Kołmogorowa.

Jak zwykle w modelu tym rozpatrujemy tylko takie przypadki że $n(x, t) \ge 0$. Należy jeszcze zadać:

1. rozkład początkowy populacji $n(x, t= 0)= n_0(x)$

2. warunki brzegowe dla $n(\pm\infty, t)$

Warunki brzegowe sformułujemy później.

Analiza równania Fishera-Kołmogorowa

1. Zauważmy, że funkcje stałe $n(x,t) = n_{0} = 0$ oraz $n(x, t) = n_{1} = 1$ są rozwiązanami tego równania. Są to te same stany stacjonarne, które są obecne w przypadku modelu Verhulsta bez dyfuzji.

2. Równanie powyższe można przeskalować w następujący sposób:

$$\tau = rt, \quad \quad y^2=\frac{r}{D} x^2, \quad \quad n(x,t) \equiv c(y, \tau)$$

Przeskalowana koncentracja spełnia równanie

$$\frac{\partial c }{\partial \tau} = c [1- c ] + \frac{\partial^2 c}{\partial y^2}, \quad \quad c=c(y, \tau)$$

Funkcje $c=0$ oraz $c=1$ są rozwiązaniami stacjonarnymi.

3. Szukamy rozwiązań w postaci fali biegnącej:

$$c(y, \tau) = U(z), \quad \mbox{gdzie } \quad z=y-v_0 \tau$$

Oznacza to, że fala przesuwa się w prawo z prędkością $v_0$. Często mówi się, że front falowy koncentracji porusza się w prawo z prędkością $v_0$.

Zauważamy, że

$$\frac{\partial c }{\partial \tau} = \frac{\partial U}{\partial \tau} = \frac{d U}{d z} \frac{\partial z}{\partial \tau} = -v_0 \frac{d U}{d z}$$

$$\frac{\partial c }{\partial y} = \frac{\partial U }{\partial y} = \frac{d U }{d z} \frac{\partial z }{\partial y} = \frac{d U }{d z}$$

$$\frac{\partial ^2 c }{\partial y^2} = \frac{\partial^2 U }{\partial y^2} = \frac{d ^2U }{d z^2}$$

Stąd, nowa funkcja $U(z)$ spełnia równanie różniczkowe zwyczajne drugiego rzędu:

$$U'' + v_0 U' + U(1-U)=0$$

Jest ono równoważne układowi 2 równań:

$$U'=V = F(U, V), \quad \quad V'=-v_0 - U(1-U) = G(U, V)$$

Stany stacjonarne określone są przez pierwiastki równań:

$$F(U, V) = 0, \quad \quad G(U, V) = 0,$$

stąd otrzymujemy 2 par rozwiązań

$$(U_0, V_0) = (0, 0), \quad (U_1, V_1) = (1, 0)$$

1. Wyznaczamy macierz Jacobiego

$$\begin{split}\quad \quad \quad \quad J = \begin{bmatrix}\frac{ \partial F}{\partial U}& \frac{\partial F}{\partial V}\\ \frac{\partial G}{\partial U}& \frac{\partial Gg}{\partial V} \end{bmatrix}= \begin{bmatrix}0, & 1\\ -1+2U, & -v_0 \end{bmatrix}\end{split}$$

w punktach stacjonarnych:

$$\begin{split}J_0= J(0, 0) = \begin{bmatrix}0& 1\\ -1& -v_0 \end{bmatrix}, \quad \quad J_1= J(1, ) = \begin{bmatrix}0& 1\\ 1& -v_0 \end{bmatrix}\end{split}$$

2. Wyznaczamy wartości własne macierzy Jacobiego $|J-\lambda I|=0$:

1. dla $(0, 0)$ otrzymujemy:

$$\lambda_{\pm}(0, 0) = (1/2)[-v_0\pm\sqrt{v_0^2-4}]$$

Należy rozważyć 2 przypadki: $v_0 \lt 2$ oraz $v_0 \ge 2$. Pierwszy przypadek należy odrzucić, ponieważ wartości własne są zespolone, o rzeczywistej części ujemnej i krzywe fazowe tworzą spiralę wokół punktu $(0, 0)$. Ale to oznacza, że $U$ przyjmuje także ujemne wartości, co nie jest dopuszczalne (koncentracje populacji mogą być tylko dodatnie lub zerowe). Pozostaje tylko drugi przypadek:

$$v_0 \ge 2$$

2. dla $(1, 0)$ otrzymujemy: $\lambda_{\pm}(1, 0) = (1/2)[-v_0\pm\sqrt{v_0^2+4}]$. Ponieważ jedna z wartości własnych jest dodatnia, ten stan stacjonarny jest niestabilny (jest to stan stacjonarny typu siodło).

sage: F(u,v)=v
sage: G(u,v)= -2*v-u+u^2 #ewolucja z  otoczenia stanu stacjonarnego (1,0)
sage: T1 = srange(0,5,0.01) #druga liczba=ilość iteracji="czas"
sage: T2 = srange(0,29,0.01)
sage: solo2=desolve_odeint(vector([F,G]), [0.8, 0.01], T1, [u,v]) # warunek pocz (U, V)
sage: solo3=desolve_odeint(vector([F,G]), [0.999, 0], T2, [u,v]) # warunek pocz (U, V)
sage: list_plot(solo2.tolist(), plotjoined=1, color='green', figsize=(7,3)) + list_plot(solo3.tolist(), plotjoined=1, figsize=(7,3))

Zauważmy, że istnieją stany początkowe z otoczenia stanu $(1, 0)$, które dążą do stanu $(0, 0)$. Innymi słowy funkcja $U(z)$ z otoczenia $U(z) \approx 1$ łączy się z funkcją $U(z)=0$. To sugeruje wybór warunków brzegowych w postaci:

$$\lim_{z\to -\infty} U(z) = 1, \quad \quad \lim_{z\to \infty} U(z) = 0$$

Pamiętając o własnościach stanów stacjonarnych w standardowym modelu Verhulsta ($n(t) \to 1$ dla $t\to \infty$), oznacza to, że dla ustalonej wartości położenia $y$:

$$\lim_{z\to -\infty} U(z) = \lim_{\tau \to \infty} U(y-v_0 \tau ) = \lim_{\tau \to \infty} c(y, \tau ) = 1$$

co jest konsystentne z zachowaniem się układu w przypadku bez migracji.

Dokładne rozwiązanie w szczególnym przypadku

Istnieje jeden szczególny przypadek, gdy znane jest dokładne rozwiązanie równania Fishera-Kołmogorowa w postaci analitycznej. Dal prędkości fali

$$v_0 = \frac{5}{\sqrt{6}}$$

rozwiązanie ma postać:

$$c(y, \tau) = \frac{1}{\left\{1+A \exp\left[ (y-v_0 \tau)/\sqrt{6}\right] \right\}^{2}}$$

gdzie stała $A\gt 0$ określa początkowy profil koncentracji.

xxxxxxxxxx

 

1

sage: A=1                                              ## określa początkowy profil koncentracji

2

sage: b=5/sqrt(6)                                      ## prędkość frontu falowego

3

sage: cs(x,t) = 1/(1+ A*exp(x-b*t)/sqrt(6))^2          ## szczególna, ale analityczna postać frontu koncentracji

4

sage: @interact

5

sage: def _(t=slider(0,5,0.2)):

6

...       pr3 = plot( cs(x,0),(x,-10,10),color='black' )

7

...       pr4 = plot( cs(x,t),(x,-10,10),fill=True )

8

...       (pr3+pr4).show(figsize=5)

Evaluate

Language: Sage Gap GP HTML Macaulay2 Maxima Octave Python R Singular

---

Messages

Analiza numeryczna równania Fishera-Kołmogorowa:

$$\frac{\partial n(x, t)}{\partial t} = r n(x, t) \;[1- n(x, t)] + D \frac{\partial^2 n(x, t)}{\partial x^2}$$

Przyjęto okresowe warunki brzegowe, poprzez odpowiedni wybór dyskretnego operatowa dyfuzji $L$. Najpierw wybieramy dyskretyzację przestrzeni, potem dobieramy maksymalny krok czasowy tak by warunek CFL był spełniony.

Aby uzyskac rozwiązanie z wyraźną falą biegnącą należy wziąć dla $r=0$ wartości dyfuzji poniżej $D\lt 0.001$. Dla większych wartości stałych dyfuzji najpierw nastąpi wyrównanie poziomu $u$ na obszarze, a potem wzrost do wartości stabilnej $u=1$.

Długość układu wchodzi efektywnie w skalowanie stałej dyfuzji, można więc bez straty ogólności rozważać układy o długości 1.

$$u^{i+1} = u^i + dt \left( r u^i (1-u^i) + D \frac{1}{h^2} Lu^i\right),$$

gdzie $L$ jest dyskretnym operatorem Laplace’a.

Ponieważ część nieliniowa jest lokalna (nie zawiera pochodnych), to warunki brzegowe są zawarte w operatorze $L$.

xxxxxxxxxx

 

1

sage: import numpy as np

2

sage: Dyf = 1.0

3

sage: r = 1.0

4

sage: l = 100.0 # dlugosc ukladu

5

sage: t_end = 50 # czas symulacje

6

sage: N = 100 # dyskretyzacja przestrzeni

7

sage: h = l/(N-1)

8

sage: dt = 0.2/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy

9

sage: sps = int(1/dt) # liczba krokow na jednostke czasu

10

sage: Nsteps=sps*t_end  # calkowita liczba krotkow

11

sage: print "sps=",sps,"dt=",dt

12

sage: one = np.identity(N)

13

sage: L=np.roll(one,-1)+np.roll(one,1)-2*one

14

sage: L[0,0]=-2.

15

sage: L[-1,-1]=-2.

16

sage: L[0,-1]=1.

17

sage: L[-1,0]=1.

18

sage: # warunek poczatkowy

19

sage: u = np.zeros(N)

20

sage: u[N/2-10:N/2+10]=.1 # small bump

21

sage: #u[:N/2]=1 # step

22

sage: Tlst=[]

23

sage: for i in range(Nsteps):

24

...       if not i%sps:

25

...           Tlst.append(list(u))

26

...       u = u + dt*(r*u*(1-u) + Dyf*(N-1)**2/l**2*L.dot(u))

27

sage: @interact

28

sage: def _(ti=slider(range(len(Tlst)))):

29

...       print r"tau=",dt*ti

30

...       p =  list_plot(Tlst[ti],plotjoined=True)

31

...       p += list_plot(Tlst[-1],plotjoined=True,color='red',ymin=0,ymax=1.5)

32

...       p += list_plot(Tlst[0],plotjoined=True,color='gray')

33

...       p.show(figsize=(8,3))

Evaluate

Language: Sage Gap GP HTML Macaulay2 Maxima Octave Python R Singular

---

Messages

xxxxxxxxxx

 

1

sage: import numpy as np

2

sage: Dyf = 1.0

3

sage: r = 1.0

4

sage: l = 100.0 # dlugosc ukladu

5

sage: t_end = 100 # czas symulacje

6

sage: N = 200 # dyskretyzacja przestrzeni

7

sage: h = l/(N-1)

8

sage: dt = 0.052/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy

9

sage: sps = int(1/dt) # liczba krokow na jednostke czasu

10

sage: Nsteps=sps*t_end  # calkowita liczba krotkow

11

sage: print "sps=",sps,"dt=",dt

12

sage: one = np.identity(N)

13

sage: L=np.roll(one,-1)+np.roll(one,1)-2*one

14

sage: L[0,0]=1.

15

sage: L[-1,-1]=1.

16

sage: # warunek poczatkowy

17

sage: u = np.zeros(N)

18

sage: #u[:int(N/2]=1 # step

19

sage: for i in range(1,3):

20

...       u[i] = 1.0 - i/3.0

21

sage: def essential_boundary_conditions(u):

22

...       u[0] = 1.

23

...       u[-1] = 0.0

24

sage: Tlst=[]

25

sage: essential_boundary_conditions(u)

26

sage: for i in range(Nsteps):

27

...       if not i%sps:

28

...           Tlst.append(list(u))

29

...       u = u + dt*(r*u*(1-u) + Dyf*(N-1)**2/l**2*L.dot(u))

30

...       essential_boundary_conditions(u)

31

sage: @interact

32

sage: def _(ti=slider(range(len(Tlst)))):

33

...       print r"t=",dt*ti

34

...       p =  list_plot(Tlst[ti],plotjoined=True)

35

...       p += list_plot(Tlst[-1],plotjoined=True,color='red',ymin=0,ymax=1.5)

36

...       p += list_plot(Tlst[0],plotjoined=True,color='gray')

37

...       p.show(figsize=(8,3))

Evaluate

Language: Sage Gap GP HTML Macaulay2 Maxima Octave Python R Singular

---

Messages

Prędkość frontu fali:

sage: pos_lst = []
sage: for T_ in Tlst:
...       for (i,a),b in zip(enumerate(T_),T_[1:]):
...           if a>=0.5 and b<=0.5:
...               pos_lst.append( i+(a-0.5)/(a-b) )
...
sage: list_plot( [l/(N-1)*(b-a)/(sps*dt) for a,b in zip(pos_lst,pos_lst[1:])] , figsize=(7,3),gridlines=[[],[2]],ymax=2)