Równanie dyfuzji w jednym wymiarze przyjmuje następującą postać:
$$ \frac{\partial u}{\partial t}= D \frac{\partial^2u}{\partial x^2}$$
Aby określić funkcje $u(x,t)$ musimy znać stan początkowy i warunku brzegowe.
$$ u(x=a,t=0)=u_a \textbf{ oraz } u(x=b,t=0)=u_b$$
lub
$$ u'(x=a,t=0)=J_a \textbf{ oraz } u'(x=b,t=0)=J_b$$
In [157]:
Out[157]:
Zaczynamy od dyskretnej postaci operatora Laplace'a:
Można przestawić go jako złożenie dwóch operatorów różnic skończonych, najlepiej w przód i w tył, tak by efekt był symetryczny.
In [2]:
N=7
B=-identity_matrix(N)
for i in range(1,B.ncols()):
B[i-1,i]=1
A=identity_matrix(N)
for i in range(1,A.ncols()):
A[i,i-1]=-1
table([[B,A,"=",B*A]])
Out[2]:
In [4]:
f = vector( [var('f%d'%i) for i in range(N)])
table([[B,f.column(),(B*f).column()]])
Out[4]:
In [6]:
table([[A,f.column(),(A*f).column()]],frame=False)
Out[6]:
In [11]:
table([[A*B,f.column(),(A*B*f).column()]])
Out[11]:
In [274]:
Out[274]:
Można też wykonać bezpośrednią konstrukcję korzystając z przybliżenia drugiej pochodnej:
$$ \displaystyle \frac{f(x_{i+1})-2 f(x_i)+f(x_{i-1}) }{h^2}\simeq \frac{\partial^2f(x)}{\partial x^2}(x_i)$$
Możemy napisać:
$$\nabla^2 f(x) = \frac{1}{h^2} L \mathbf{f}, $$
gdzie:
$$ L = \left(\begin{array}{rrrr} -2.0 & 1.0 & 0.0 & 0.0 \\ 1.0 & -2.0 & 1.0 & 0.0 \\ 0.0 & 1.0 & -2.0 & 1.0 \\ 0.0 & 0.0 & 1.0 & -2.0 \end{array}\right)$$
a $\mathbf{f}$ jest wektorem wartości funkcji wypróbkowanym na dziedzinie.
In [12]:
N=7
L = matrix(RDF,N)
for i in range(1,N-1):
L[i,i-1],L[i,i],L[i,i+1] = 1,-2, 1
L[0,0],L[0,1],L[0,-1] = -2, 1, 1
L[-1,-1],L[-1,-2],L[-1,0] = -2, 1, 1
show(L[:10,:10])
Można też skorzystać z definicji i obliczyć przybliżenie pochodnej stosując technikę "array slicing". Na przykład używając wydajnych numerycznie tablic numpy mamy:
In [14]:
import numpy as np
u = np.random.randint(5,size=10)*1.0
a = np.zeros_like(u)
b = np.zeros_like(u)
c = np.zeros_like(u)
d = np.zeros_like(u)
a[1:-1]=u[2:]-2.0*u[1:-1]+u[:-2]
b[1:-1]=u[2:]
c[1:-1]=-2.0*u[1:-1]
d[1:-1]=u[:-2]
table([[u],[d],[b],[c],'=',[a]])
Out[14]:
Dalej, będziemy dla wygody stosować skonczony operator liniowy $L$. Zobaczny jak działa taki operator na dolowlną funkcję wypróbkowaną na równoodległych punktach
In [16]:
f = vector( [var('f%d'%i) for i in range(N)])
table([[L,'$\cdot$',f.column() ,'$=$', (L*f).column() ]])
Out[16]:
In [176]:
#print latex(L), latex(f.column()), latex((L*f).column() )
#var('dt,D')
#f = vector( [var('u_i%d'%i) for i in range(N)])
#print latex(identity_matrix(N)),latex(dt*D), latex(L), latex(f.column()), latex((L*f).column() )
Out[176]:
In [178]:
Out[178]:
Równanie dyfuzji możemy zdyskretyzować w dziedzinie czasowej stosując jawny lub niejawny schemat Eulera i zastępując Laplace'an jego dyskretnym odpowiednikiem:
$$u^{i+1} = u^i + \left[ \frac{D dt}{h^2} \right] Lu^i.$$
Schemat, jak i jego zbieżność zależą od jednej stałej $\frac{D dt}{h^2}$, zwanej też liczbą Couranta, w której tkwią zarówno wielkości fizyczne jak i dyskretyzacja układu. Dalej, zapisując układ iteracji jako:
$$u^{i+1} = \left( I + L\frac{D dt}{h^2} \right) u^i,$$
widzimy, że rozwiązanie jest potęgami $L_t^1,L_t^2,L_t^3,\dots,L_t^n$ operatora
$$L_t= I + \frac{D dt}{h^2} L.$$
Schemat niejawny zawiera wartość $u$ w chwili $i+1$ po prawej stronie:
$$u^{i+1} = u^i + \left[ \frac{D dt}{h^2} \right] Lu^{i+1},$$
co przepisując szukaną $u_{i+1}$ na prawą stronę daje nam niejednorodny układ rówń liniowych:
$$ \left( I - \frac{D dt}{h^2} L \right) u^{i+1}= u^{i}.$$
W przypadku małego $dt$ schematy te są równoważne. Korzystając z twierdzenia o macierzowym szeregu geometrycznym, lub rozwinięcia w szereg Taylora funkcji macierzowej (link) mamy:
$$ \left({I -\epsilon A}\right)^{-1} = I+\epsilon A+ \epsilon^2 A^2 \dots $$
Widzimy, że formalnie rozwiązując układ równań liniowych ze schematu niejawnego, biorąc wyrazy liniowe w $dt$ dostajemy schemat jawny.
In [17]:
u = vector(RDF,[0,0,0,1,0,0,0])
print u
print L*u
In [228]:
Out[228]:
Dyskretny operator Laplace'a w 1d potrzebuje wartości funkcji z dwóch węzłów sąsiadujących z węzłem w którym obliczamy wartosc laplasjanu. Jeśli punkt ten znajduje się na brzegu obszaru to brakuje informacji by obliczyć poprawny operator. Informację tą trzeba dostarczyc w postaci warunku brzegowego by problem miał jednoznaczne rozwiązanie. Przykładem warunku brzegowego może być:
In [183]:
Out[183]:
Załóżmy, że nakładamy warunek Dirichleta na końcach przedziału: $u_0=1$ i $u_{N-1}=2$. Oznacza to, że w ewolucji czasonej wektora $u$, wartości na brzegu będą zawsze utrzymywane "zewnetrznie" na zadanych wartościach. Biorąc jawny schemat mamy:
$$
\left(\begin{array}{r}\mathbf{u_0^{i+1}}\\u_1^{i+1}\\u_2^{i+1}\\u_3^{i+1}\\\mathbf{u_4^{i+1}}\end{array}\right) =
\underbrace{
\left[ \left(\begin{array}{rrrrr}\mathbf{ 1 }& 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array}\right) +\frac{dt D}{h^2} \left(\begin{array}{rrrrr} -2.0 & 1.0 & 0.0 & 0.0 & 0.0 \\ 1.0 & -2.0 & 1.0 & 0.0 & 0.0 \\ 0.0 & 1.0 & -2.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 1.0 & -2.0 & 1.0 \\ 0.0 & 0.0 & 0.0 & 1.0 & -2.0 \end{array}\right)\right]
}_{L_t}
\left(\begin{array}{r}\mathbf{u_0^i}\\u_1^i\\u_2^i\\u_3^i\\\mathbf{u_4^i}\end{array}\right)
$$
otrzymamy po każdym kroku wartości $u_0^{i+1}$ i $u_4^{i+1}$, które nie będą spełniały dokładnie warunku brzegowego. Należy wieć po każdym kroku wymusić wartości:
$$
u_0^{i+1}=1 \quad u_4^{i+1}=2
$$
Zauważmy, że wtedy w kolejnym kroku wartościami krańcowymi wektora $u$ będą wielkości zgodne z warunkami i będą dawały poprawny przyczynek do swoich sąsiadów.
Warunki Dirichleta, są zwane "istotnymi warunkami brzegowymi" (essential boundary conditions). Nie da się tak zmodyfikować dyskretnego operatora $L$ by automatycznie spełniał te warunki i muszą być dołączone dodatkowo w schemacie numerycznym.
In [18]:
# Dirichlet
def init_L(N=7):
L = matrix(RDF,N)
for i in range(1,N-1):
L[i,i-1],L[i,i],L[i,i+1] = 1,-2, 1
L[0,0],L[-1,-1] = 1, 1
return L
def essential_boundary_conditions(u):
u[0] = 1.2
u[-1] = 2.1
L = init_L(7)
show(L[:10,:10])
Okresowy warunek brzegowy w przypadku jednowymiarowym polega na utożsamieniu $u_0=u_{N}$. Obszar na którym rozwiązywane jest równanie jest topologicznie równoważny okręgowi. Okrąg nie posiada brzegu więc problem jest dobrze określony - nie ma gdzie zadawać warunku brzegowego.
Warunek ten można zaimplementowac modyfikująć dyskretny operator Laplace'a $L$ tak by:
$L_{0,N-1}=1$ i $L_{N-1,0}=1$. Niech $N=5$, mamy:
$$ \left(\begin{array}{rrrrr} -2.0 & 1.0 & 0.0 & 0.0 & 1.0 \\ 1.0 & -2.0 & 1.0 & 0.0 & 0.0 \\ 0.0 & 1.0 & -2.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 1.0 & -2.0 & 1.0 \\ 1.0 & 0.0 & 0.0 & 1.0 & -2.0 \end{array}\right) \left(\begin{array}{r} f_{0} \\ f_{1} \\ f_{2} \\ f_{3} \\ f_{4} \end{array}\right) \left(\begin{array}{r} -2.0 \, f_{0} + f_{1} + f_{4} \\ f_{0} - 2.0 \, f_{1} + f_{2} \\ f_{1} - 2.0 \, f_{2} + f_{3} \\ f_{2} - 2.0 \, f_{3} + f_{4} \\ f_{0} + f_{3} - 2.0 \, f_{4} \end{array}\right) $$
Widać, że taki operator oblicza poprawnie Laplacjan dla punktów skrajcym, biarąc za brakujące punkty $u_{-1}$ i $u_5$, odpowiednio: $u_{4}$ oraz $u_0$.
In [19]:
# PBC
def init_L_pbc(N=7):
L = matrix(RDF,N)
for i in range(1,N-1):
L[i,i-1],L[i,i],L[i,i+1] = 1,-2, 1
L[0,0],L[0,1],L[0,-1] = -2, 1, 1
L[-1,-1],L[-1,-2],L[-1,0] = -2, 1, 1
return L
def essential_boundary_conditions(u):
pass
L = init_L_pbc(7)
show(L[:10,:10])
In [225]:
L.rank()
Out[225]:
In [231]:
Out[231]:
W przypadku ogólnym, rozważmy równania dające się zapisać w postaci prawa zachowania:
$$\frac{\partial u}{\partial t} = - \nabla \cdot \vec J,$$
gdzie $J$ to strumień pola $u$. Równanie dyfuzji można przedstawić z tej postaci przy założeniu że:
$$\vec J = - \vec\nabla u$$
Jeśli równanie zawiera człon adwekcyjny (tzn proporcjonalny do pierwszej pochodnej) to strumień będzie zawierał dodatkowe człony.
Widać, że przypadku jednowymiarowego równania dyfizji warunek von Neumanna jest efektywnie warunkiem na pochodną funkcji na brzegu:
$$\frac{u_1-u_0}{h}=-J$$
Przypadkiem szczególnym warunku Neumanna jest bariera odbijająca, w której zakładamy że strumień cząstek opisywanych gęstością lub stężeniem $u$ przez barierę wynosi zero. W takim przypadku można napisać operator Laplace'a, który będzie konsystentny z tym warunkiem:
$$
\left(\begin{array}{rrrrr} -1.0 & 1.0 & 0.0 & 0.0 & 0.0 \\ 1.0 & -2.0 & 1.0 & 0.0 & 0.0 \\ 0.0 & 1.0 & -2.0 & 1.0 & 0.0 \\ 0.0 & 0.0 & 1.0 & -2.0 & 1.0 \\ 0.0 & 0.0 & 0.0 & 1.0 & -1.0 \end{array}\right) \left(\begin{array}{r} f_{0} \\ f_{1} \\ f_{2} \\ f_{3} \\ f_{4} \end{array}\right)
= \left(\begin{array}{r} -f_{0} + f_{1} \\ f_{0} - 2.0 \, f_{1} + f_{2} \\ f_{1} - 2.0 \, f_{2} + f_{3} \\ f_{2} - 2.0 \, f_{3} + f_{4} \\ f_{3} - f_{4} \end{array}\right) $$
Widać, że taki operator, zamiast drugiej pochodnej w punktach skrajnych oblicza pierwszą pochodną. Intuicyjnie, działanie operatora ewolucji na dowolny wektor będzie poprawiało wartość w pierwszym i ostatnim węźle tak długo aż pierwsze pochodne będą zero.
Warto odnotować, że taki operator ma rząd o jednej mniejszy od wymiaru. Wynika z tego, że rozwiązanie zerowe spełnia takie równanie. Rzeczywiście: równanie dyfuzji na obszarze z odbijającymi scianami jest spełnione jeśli w układzie nie ma cząstek! Ponadto widać, że jesli rozwiązanie jest określone co do wartości stałej multyplikatywnej.
In [20]:
# von Neumann/reflecting BC
def init_L_ref(N=7):
L = matrix(RDF,N)
for i in range(1,N-1):
L[i,i-1],L[i,i],L[i,i+1] = 1,-2, 1
L[0,0],L[0,1] = -1, 1
L[-1,-1],L[-1,-2] = -1, 1
return L
def essential_boundary_conditions(u):
pass
L = init_L_ref(7)
show(L[:10,:10])
In [226]:
L.rank()
Out[226]:
In [168]:
Out[168]:
In [21]:
Lt=identity_matrix(N)+0.40*L
eig = list(Lt.eigenvalues())
eig_s = sorted(map(lambda x:x.n(digits=3),map(real,eig)))
show(eig_s)
Sprawdźmy wartości własne operatora $L_t=I+\frac{D dt}{h^2}L$, dla różnych wartości stałej $C=\frac{D dt}{h^2}$. Zacznijmy od małej wartości np: $C=0.2$. Dla $N=5$ i operatora z okresowymi warunkami brzegowymi otrzymujemy:
$$\left[0.240, 0.240, 0.511, 0.511, 0.849, 0.849, 1.00\right].$$
Widać, że wartości własne są rzeczywiste, dodatnie, mniejsze od jednego z wyjątkiem jednej. Ewolucja czasowa układu jest dana przez potęgi operatora $L_t$: $$L_t^1,L_t^2,L_t^3,\dots,L_t^n.$$ Oznacza to, że kolejne iteracje będą wygaszać składowe wektora wzdłuż wszystkich wektorów własnych, z wyjątkiem tego należącego do wartości jeden, która to będzie stanem stacjonarnym.
Niech $C=0.4$, otrzymujemy wtedy:
$$\left[-0.521, -0.521, 0.0220, 0.0220, 0.699, 0.699, 1.00\right]$$
Pojawiają się ujemne wartości własnych, co oznacza oscylacje pomiędzy dodatnimi i ujemnymi wartościami np. stężenia w czasie. Nie jest to efekt fizyczny i jawny algorytm Eulera dla równaia dyfuzji dla $C=0.4$ jest robieżny.
Warto odnotować, że stała od której zależy stabilnośc zawiera w liczniku krok czasowu a w mianowniku kwadrat kroku przestrzennego. Oznacza to, że zmniejszając dyskretyzacje przestrzeni musimy jednocześnie używać mniejszego kroku czasowego, by schemat był stabilny.
In [22]:
N=7
L = matrix(RDF,N)
for i in range(1,N-1):
L[i,i-1],L[i,i],L[i,i+1] = 1,-2, 1
L[0,0],L[0,1],L[0,-1] = -2, 1, 1
L[-1,-1],L[-1,-2],L[-1,0] = -2, 1, 1
@interact
def _(C=slider(0.01,1.0,0.01)):
Lt=matrix(RDF,identity_matrix(N)+C*L)
eig = list(Lt.eigenvalues())
l = sorted(map(lambda x:x.n(digits=3),map(real,eig)))
print l[0:6],"...",l[-1]
In [182]:
Out[182]:
Mając juz wszystkie składniki można napisać algorytm który będzie rozwiązywał numerycznie równanie dyfuzji przy zadanych warunkach brzegowych i początkowych.
In [181]:
Out[181]:
In [180]:
L.ncols(),L.rank()
Out[180]:
In [232]:
Out[232]:
In [162]:
L = init_L_ref(45)
def essential_boundary_conditions(u):
pass
Tlst=[]
Lt=matrix(RDF,identity_matrix(L.ncols())+0.2*L)
u = zero_vector(RDF,L.ncols())
u[ int(L.ncols()/2) ] = 1.0
essential_boundary_conditions(u)
for i in range(150):
Tlst.append(u)
u = Lt*u # schemat jawny
essential_boundary_conditions(u)
@interact
def _(ti=slider(range(len(Tlst)))):
p = list_plot(Tlst[ti],plotjoined=True)
p += list_plot(Tlst[-1],plotjoined=True,color='gray',ymin=-0.2,ymax=1.0)
p += list_plot(Tlst[0],plotjoined=True,color='gray')
p.show(figsize=(9,3))
Out[162]:
Warunek unormowania:
In [161]:
[sum(T_) for T_ in Tlst]
Out[161]:
Rozważmy równanie:
$$ \frac{\partial u}{\partial t}= D \frac{\partial^2u}{\partial x^2}$$
na odcinku $(0,l)$ z odbijającymi warunkami brzegowymi. W tym celu stosujemy jawny schemat Eulera. Krok przestrzenny $h$ jest równy:
$$h = \frac{l^2}{(N-1)^2}.$$
Wobec tego mamy następujący infinitezymalny operator ewolucji
$$L_t= I + dt\frac{D (N-1)^2}{l^2} L u^i,$$
przy czym maksymalny krok czasowy zależy od parametrów układu i jest ograniczony przez:
$$dt_{max}<0.25 \frac{l^2}{(N-1)^2 D}. $$
In [259]:
N = 55
Dyf = 1.0
l = 100.
dx = l/(N-1)
dt_max = 0.2*dx^2/Dyf
dt = dt_max/2.0
C = dt*Dyf/dx^2
Nsteps = 1200
print C,dt,dx
Out[259]:
In [276]:
print "Czas dyfuzyjny:",l^2/Dyf/dt
Out[276]:
In [277]:
Nsteps = int(l^2/Dyf/dt)
Out[277]:
In [154]:
L = init_L_ref(N)
def essential_boundary_conditions(u):
pass
Tlst=[]
Lt=matrix(RDF,identity_matrix(L.ncols())+C*L)
u = zero_vector(RDF,L.ncols())
u[ int(L.ncols()/5) ] = 1.0/dx
essential_boundary_conditions(u)
for i in range(Nsteps):
Tlst.append(u)
u = Lt*u # schemat jawny
essential_boundary_conditions(u)
@interact
def _(ti=slider(range(len(Tlst)))):
p = list_plot(Tlst[ti],plotjoined=True)
p += list_plot(Tlst[-1],plotjoined=True,color='gray',ymin=-0.2,ymax=0.5)
p += list_plot(Tlst[0],plotjoined=True,color='gray')
p.show(figsize=(9,3))
Out[154]:
In [261]:
Out[261]:
In [262]:
c(x,t)=1/sqrt(4*pi*Dyf*t)*exp(-(x^2)/(4*Dyf*t) )
print "Unormowanie wzoru analitycznego:",integrate(c(x,0.23),(x,-oo,oo))
T = [i*dt for i in range(150)]
X = [ (-l/5. + i*dx).n() for i in range(N)]
@interact
def _(ti=slider(range(1,len(Tlst)))):
print "t=",dt*ti,"Norma=",sum(Tlst[ti])*(l/(N-1))
plt = point(zip(X,Tlst[ti]),figsize=(7,3),color='red')
plt += plot(c(x,dt*ti),(x,-l/2,l/2))
plt.show(figsize=(8,3))
Out[262]:
In [132]:
Out[132]:
In [70]:
Out[70]:
In [210]:
import numpy as np
from scipy.sparse import dia_matrix
Out[210]:
In [211]:
%timeit
sparse = True
slicing = False
Dyf = 1.0
r = 1.0
l = 100.0 # dlugosc ukladu
t_end = 100 # czas symulacje
N = 250 # dyskretyzacja przestrzeni
h = l/(N-1)
dt = 0.052/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy
sps = int(1/dt) # liczba krokow na jednostke czasu
Nsteps=sps*t_end # calkowita liczba krotkow
print "sps=",sps,"dt=",dt,'Nsteps=',Nsteps
if sparse:
L = dia_matrix( (np.array([N*[-2.],N*[1.],N*[1.]]),np.array([0,-1,1])), shape=(N,N))
if slicing:
one = np.identity(N)
L=np.roll(one,-1)+np.roll(one,1)-2*one
L[0,0]=1.
L[-1,-1]=1.
# warunek poczatkowy
u = np.zeros(N)
#u[int(N/2)-20:int(N/2)+20]=1 # step
#for i in range(1,3):
# u[i] = 1.0 - i/3.0
u[:int(N/2)]=1
def essential_boundary_conditions(u):
u[0] = 1.0
u[-1] = 0.0
Tlst=[]
essential_boundary_conditions(u)
for i in range(Nsteps):
if not i%sps:
Tlst.append(list(u))
if slicing:
u[1:-1] = u[1:-1] + dt*(r*u[1:-1]*(1-u[1:-1]) + Dyf*(N-1)**2/l**2*np.diff(u,2))
else:
u = u + dt*(r*u*(1-u) + Dyf*(N-1)**2/l**2*L.dot(u))
essential_boundary_conditions(u)
print "Saved ",len(Tlst), " from ", Nsteps
Out[211]:
In [266]:
pos_lst = []
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) )
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)
Out[266]:
In [235]:
@interact
def _(ti=slider(range(len(Tlst)))):
print r"t=",dt*ti
p = list_plot(Tlst[ti],plotjoined=True)
p += list_plot(Tlst[-1],plotjoined=True,color='red',ymin=0,ymax=1.5)
p += list_plot(Tlst[0],plotjoined=True,color='gray')
p.show(figsize=(8,3))
Out[235]:
In [237]:
%timeit
Dyf = 1.0
r = 1.0
l = 100.0 # dlugosc ukladu
t_end = 100 # czas symulacje
N = 100 # dyskretyzacja przestrzeni
h = l/(N-1)
dt = 0.052/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy
sps = int(1/dt) # liczba krokow na jednostke czasu
Nsteps=sps*t_end # calkowita liczba krotkow
print "sps=",sps,"dt=",dt,'Nsteps=',Nsteps
# warunek poczatkowy
u = np.zeros((N,N))
#u[int(N/2)-5:int(N/2)+5,int(N/2)-5:int(N/2)+5]=1 # step
#u[:int(N/2)+5,:]=1.0 # step
u[int(N/2),int(N/2)]=1.0
def essential_boundary_conditions(u):
u[:,0] = 0.0
u[:,-1] = 0.0
u[-1,:] = 0.0
u[0,:] = 0.0
Tlst=[]
essential_boundary_conditions(u)
for i in range(Nsteps):
if not i%sps:
Tlst.append(u.copy())
u[1:-1,1:-1] = u[1:-1,1:-1] + dt*(r*u[1:-1,1:-1]*(1-u[1:-1,1:-1]) + \
Dyf*(N-1)**2/l**2*(np.diff(u,2,axis=0)[:,1:-1]+np.diff(u,2,axis=1)[1:-1,:]))
essential_boundary_conditions(u)
print "Saved ",len(Tlst), " from ", Nsteps
Out[237]:
In [215]:
import pylab
@interact
def _(ti=slider(range(len(Tlst)))):
print r"t=",dt*ti
if True:
pylab.clf()
pylab.imshow(Tlst[ti],origin='top')
pylab.savefig('1.png',dpi=70)
else:
p = matrix_plot(Tlst[ti])
p.show(figsize=(4,4))
Out[215]:
In [268]:
Out[268]:
Dynamika modelu bez dyfuzji.
In [267]:
a=1.0
b=0.1
eps=0.1
a = 0.75
b = 0.0006
eps = 0.072
var('u v')
f(u,v) = u*(1-u)*(u-(v-b)/a)
g(u,v) = u-v
V = vector( (1/eps*f,g))
V=V/V.norm()
vfield=plot_vector_field(V,(u,0,1),(v,0,1))+implicit_plot(g,(u,0,1),(v,0,1))
t = srange(0,4/eps,0.01)
sol = desolve_odeint([19*f,g], [0.5,0.0], t, [u,v])
plt_phase = vfield+line(sol,color='red',figsize=5)
plt_time = line(zip(t,sol[:,0]),figsize=5)
html.table([[plt_phase,plt_time]])
Out[267]:
In [270]:
Out[270]:
In [218]:
%timeit
import numpy as np
sparse = True
slicing = True
Dyf_u = 1.0
Dyf_v = 0.052
Dyf = max(Dyf_u,Dyf_v)
a = 1.0
b = 0.001
eps = 0.072
l = 100.0 # dlugosc ukladu
t_end = 100 # czas symulacje
N = 160 # dyskretyzacja przestrzeni
h = l/(N-1)
dt = 0.052/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy
dt_dyn = (1.0/eps)/125.0
sps = int(1/dt) # liczba krokow na jednostke czasu
Nsteps=sps*t_end # calkowita liczba krotkow
print "dt,dt_dyn",dt,dt_dyn
dt = min(dt,dt_dyn)
print "sps=",sps,"dt=",dt,'Nsteps=',Nsteps
# warunek poczatkowy
u = np.zeros((N,N))
v = np.zeros((N,N))
#u[int(N/2)-5:int(N/2)+5,int(N/2)-5:int(N/2)+5]=1 # step
#u[:int(N/2)+5,:]=1.0 # step
#u[int(N/2)-5:int(N/2)+5,int(N/2)-20:int(N/2)+20]=1.0
#v[int(N/2)-5:int(N/2)+3,int(N/2)-20:int(N/2)+20]=1.0
#u[int(N/2)-5:int(N/2)+5,int(N/2)-20:int(N/2)+20]=1.0
#v[int(N/2)-5:int(N/2)+5,int(N/2)-22:int(N/2)+18]=1.0
#u[-20:-1,int(N/2)-5:int(N/2)+5]=1.0
#v[-20:-1,int(N/2)-6:int(N/2)+4]=1.0
#u[:5,:]=1.0
#v[:4,:]=1.0
#u[-10:,:]=1.0
#v[-4:,:]=1.0
u[:int(N/2),int(N/2)-5:int(N/2)+5]=1.0
v[:int(N/2),int(N/2)-6:int(N/2)+4]=1.0
# aby wymusic ruch falowy, przesuwamy u wzgledem v
def essential_boundary_conditions(u):
u[:,0] = 0.0
u[:,-1] = 0.0
u[-1,:] = 0.0
u[0,:] = 0.0
v[:,0] = 0.0
v[:,-1] = 0.0
v[-1,:] = 0.0
v[0,:] = 0.0
Tlst=[]
Tvlst=[]
essential_boundary_conditions(u)
for i in range(Nsteps):
if not i%sps:
Tlst.append(u.copy())
Tvlst.append(v.copy())
u[1:-1,1:-1] = u[1:-1,1:-1] + dt*(1.0/eps*u[1:-1,1:-1]*(1-u[1:-1,1:-1])*( u[1:-1,1:-1]-(v[1:-1,1:-1]+b)/a ) + \
Dyf_u*(N-1)**2/l**2*(np.diff(u,2,axis=0)[:,1:-1]+np.diff(u,2,axis=1)[1:-1,:]))
v[1:-1,1:-1] = v[1:-1,1:-1] + dt*( (u[1:-1,1:-1]-v[1:-1,1:-1]) )
# + \
# Dyf*(N-1)**2/l**2*(np.diff(v,2,axis=0)[:,1:-1]+np.diff(v,2,axis=1)[1:-1,:]))
essential_boundary_conditions(u)
print "Saved ",len(Tlst), " from ", Nsteps
Out[218]:
In [250]:
anim=animate([matrix_plot(u,cmap='jet',figsize=(2,2)) for u in Tlst[:35]])
anim.show()
Out[250]:
In [221]:
import pylab
@interact
def _(ti=slider(range(len(Tlst)))):
print r"t=",dt*ti*sps
if True:
pylab.subplot(1,2,1)
pylab.imshow(Tlst[ti],vmin=0,vmax=1,origin='top')
pylab.subplot(1,2,2)
pylab.imshow(Tvlst[ti],vmin=0,vmax=1,origin='top')
pylab.savefig('1.png',dpi=70)
else:
p = matrix_plot(Tlst[ti])
p.show(figsize=(4,4))
Out[221]:
In [244]:
anim=animate([matrix_plot(u,cmap='jet',figsize=(2,2)) for u in Tlst[40:72:2]])
anim.show()
Out[244]:
In [247]:
Out[247]:
In [271]:
Out[271]:
In [245]:
Out[245]:
In [248]:
Out[248]:
In [242]:
Out[242]:
Dynamika modelu (bez dyfuzji)
In [251]:
F = 0.1
k = 0.04
F = 0.04
k = 0.062
var('u v')
f(u,v) = -u*v^2+ F*(1-u)
g(u,v) = u*v^2-(F+k)*v
V = vector( (f,g))
V=V/V.norm()
vfield=plot_vector_field(V,(u,-0.1,1.2),(v,-0.1,1.2))
nullclines=implicit_plot(f,(u,-0.1,1.2), (v,-0.1,1.2),color='red')+implicit_plot(g,(u,-0.1,1.2),(v,-0.1,1.2),color='green')
t = srange(0,226,0.01)
plt_phase = vfield+nullclines
sol = desolve_odeint([f,g], [1,0.25], t, [u,v])
plt_phase += line(sol,color='blue',figsize=5)
plt_time = line(zip(t,sol[:,0]),figsize=5) + line(zip(t,sol[:,1]),color='green')
html.table([[plt_phase,plt_time]])
Out[251]:
In [258]:
Out[258]:
In [243]:
%timeit
Dyf_u = 2e-5
Dyf_v = 2e-5
Dyf = max(Dyf_u,Dyf_v)
F = 0.04
k = 0.0562
l = 2.5 # dlugosc ukladu
t_end = 100 # czas symulacje
N = 256 # dyskretyzacja przestrzeni
h = l/(N-1)
dt = 0.052/(Dyf*(N-1)**2/l**2) # 0.2 z warunku CFL, krok nie moze byc wiekszy
#dt_dyn = (1.0)/125.0
#dt = min(dt,dt_dyn)
#sps = int(1/dt) # liczba krokow na jednostke czasu
#Nsteps=sps*t_end # calkowita liczba krotkow
n=1
dt = 1.0/n
Nsteps = 20000*n
sps = n*10
print "dt,dt_dyn",dt,dt_dyn
print "sps=",sps,"dt=",dt,'Nsteps=',Nsteps
# warunek poczatkowy
u = np.ones((N,N))
v = np.zeros((N,N))
u[int(N/2)-10:int(N/2)+20,int(N/2)-10:int(N/2)+10]=.5
v[int(N/2)-10:int(N/2)+10,int(N/2)-20:int(N/2)+10]=.24
u = u + np.random.uniform(low=0,high=0.01,size=(N,N))
v = v + np.random.uniform(low=0,high=0.01,size=(N,N))
def essential_boundary_conditions(u):
u[:,0] = 1.0
u[:,-1] = 1.0
u[-1,:] = 1.0
u[0,:] = 1.0
v[:,0] = 0.0
v[:,-1] = 0.0
v[-1,:] = 0.0
v[0,:] = 0.0
Tlst=[]
Tvlst=[]
essential_boundary_conditions(u)
for i in range(Nsteps):
if not i%sps:
Tlst.append(u.copy())
Tvlst.append(v.copy())
u[1:-1,1:-1] = u[1:-1,1:-1] + dt*(-u[1:-1,1:-1]*v[1:-1,1:-1]**2 + F*(1.0r-u[1:-1,1:-1]) + \
Dyf_u*(N-1)**2/l**2*(np.diff(u,2,axis=0)[:,1:-1]+np.diff(u,2,axis=1)[1:-1,:]))
v[1:-1,1:-1] = v[1:-1,1:-1] + dt*( u[1:-1,1:-1]*v[1:-1,1:-1]**2 - (F+k)*v[1:-1,1:-1]+ \
Dyf_v*(N-1)**2/l**2*(np.diff(v,2,axis=0)[:,1:-1]+np.diff(v,2,axis=1)[1:-1,:]))
essential_boundary_conditions(u)
print "Saved ",len(Tlst), " from ", Nsteps
pylab.clf()
pylab.imshow(u,vmin=0.2,vmax=1,origin='top')
pylab.colorbar()
pylab.savefig('1.png',dpi=70)
Out[243]:
In [257]:
import pylab
@interact
def _(ti=slider(range(len(Tlst)))):
print r"t=",dt*ti*sps
pylab.clf()
pylab.subplot(1,2,1)
pylab.imshow(Tlst[ti],vmin=0,vmax=1,origin='top')
pylab.subplot(1,2,2)
pylab.imshow(Tvlst[ti],vmin=0,vmax=1,origin='top')
pylab.colorbar()
pylab.savefig('1.png',dpi=70)
Out[257]:
In [256]:
Out[256]:
In [255]:
Out[255]:
In [254]:
Out[254]:
In [216]:
Out[216]:
In [214]:
Out[214]:
In [208]:
Out[208]:
In [212]:
Out[212]:
In [201]:
Out[201]:
In [203]:
Out[203]:
In [200]:
Out[200]:
In [202]:
Out[202]: