# -*- coding: utf-8 -*-
import matplotlib
#matplotlib.use('TkAgg')
from scipy.integrate import odeint
from scipy.optimize import fsolve
from pylab import figure, plot, show, xlabel, ylabel ,draw, ion,ylim,ioff, pause
import numpy as np

def BZ(y,t,n):
    '''
    równania :
    n - ilość węzłów
    
    du/dt = d2u/dx2 +u*(1-u-r*v)
    dv/dt = d2v/dx2 -b*u*v 
    
    warunek początkowy: u(x,0)= 0
                        v(x,0)=0
    warunki brzegowe: 
        z lewej: 
            du/dx(0,t)  = 0
            v(0,t) = 1
        z prawej:
            u(1,t) = 1
            dv/dx(1,t) = 0
    '''
    u,v=np.hsplit(y,2)
    b=1
    r=50
    D = 1e-4
    # bierzemy siatkę o n węzłach
    h = 1./n # to jest odstęp między węzłami
    h2= h**2
    du = np.zeros(n)
    A= np.zeros((n,n))
    for wiersz in range(n):
        for kol in range(n):
            if wiersz==kol:
                A[wiersz,kol] =-2
            elif kol+1==wiersz:
                A[wiersz,kol] =1
            elif kol-1==wiersz:
                A[wiersz,kol] = 1
    A[0,0]+=1           
    A *= D/h2
    
    f = u*(1-u-r*v)
    a =np.zeros(n)
    a[-1] = 1
    a *= D/h2
    du = np.dot(A,u) + a #+f
    
    dv = np.zeros(n)
    B= np.zeros((n,n))
    for wiersz in range(n):
        for kol in range(n):
            if wiersz==kol:
                B[wiersz,kol] =-2
            elif kol+1==wiersz:
                B[wiersz,kol] =1
            elif kol-1==wiersz:
                B[wiersz,kol] = 1
    B[-1,-1]+=1           
    B *= D/h2
    
    g = -b*u*v
    wb =np.zeros(n)
    wb[0] = 1
    wb *= D/h2
    dv = np.dot(B,v) + wb #+g
    
    dy=np.hstack((du,dv))
    
    return dy




## całkowanie


# prametry

time = np.linspace(0.0,50.0,1000) 

n=100#ilość węzłów
x = np.linspace(0,1,n)
# warunki początkowe
uinit = np.zeros(n)
uinit[-1] = 1
vinit = np.zeros(n)
vinit[0] = 1
yinit = np.hstack((uinit[1:-1],vinit[1:-1]))



y = odeint(BZ, yinit, time, args =(n-2,))

u,v=np.hsplit(y,2)


ion()
lu,= plot(x[1:-1],u[0,:])
lv,= plot(x[1:-1],v[0,:])

ylim([-1,1])
for t in range(0, y.shape[0]):
    lu.set_ydata(u[t,:]) # y[:,0] pierwsza kolumna y - zawiera zmienną v
    lv.set_ydata(v[t,:])
    draw()
    pause(0.01)


xlabel('x') 
ylabel('u') 



show()