TI/Programowanie dla Fizyków Medycznych/RRZ
Z Brain-wiki
import numpy as np
import pylab as py
#rozwiazujemy rownanie dx(t)/dt=f(t,x)
#metoda Explicit Euler
#f - funkcja f z rownania
#x0-wartosc poczatkowa
#t0-czas poczatkowy
#tk-czas koncowy
#h-krok czasowy
def EE(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
#wektor wynikowy
if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
else: x=np.zeros(N)
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
#index
i=1
#petla glowna
while (i<N):
x[i]=np.array(x[i-1]+h*f(t[i-1],x[i-1]))
i+=1
return t,x
def leapfrog(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
#wektor wynikowy
if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
else: x=np.zeros(N)
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
#index
i=1
#petla glowna
while (i<N):
k1=f(t[i-1],x[i-1])
k2=f(t[i-1]+h*0.5,x[i-1]+0.5*h*k1)
x[i]=np.array(x[i-1]+h*k2)
i+=1
return t,x
def Heun(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
#wektor wynikowy
if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
else: x=np.zeros(N)
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
#index
i=1
#petla glowna
while (i<N):
k1=f(t[i-1],x[i-1])
k2=f(t[i-1]+h*0.5,x[i-1]+0.5*h*k1)
x[i]=np.array(x[i-1]+h*0.5*(k1+k2))
i+=1
return t,x
def RK4(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
#wektor wynikowy
if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
else: x=np.zeros(N)
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
#index
i=1
#petla glowna
while (i<N):
k1=h*f(t[i-1],x[i-1])
k2=h*f(t[i-1]+h*0.5,x[i-1]+0.5*k1)
k3=h*f(t[i-1]+h*0.5,x[i-1]+0.5*k2)
k4=h*f(t[i-1]+h,x[i-1]+k3)
x[i]=np.array(x[i-1]+(k1+2.0*k2+2.0*k3+k4)/6)
i+=1
return t,x
#Metoda Verlet sluzy do rozw. rownan rowniczkowych typu:
# d^2 x(t)/dt^2=a(x(t))
# x_(n+1)=x_n+v_n*h+1/2 a(x_n)*h^2
# v_n=(x_n-x_(n-1))/h
#podstawiamy do poprzedniego
# x_(n+1)=2*x_n-x_(n-1)+1/2 a(x_n)*h^2
#przekształcają otrzymamy
# x_(n-1)=2*x_n-x_(n+1)+1/2 a(x_n)*h^2
#zatem algorytm jest symetryczny w czasie
#rownanie ma postac:
#d(x,x')/dt=(x',x'')=f(t,(x,x'))
def Verlet(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
if len(x0)!=2:
print 'algorytm dziala tylko dla rown. 2 rzedu'
return None
x=np.zeros((N,2))
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
x[1]=np.array(x0)
#index
i=2
#petla glowna
while (i<N):
x_i=2*x[i-1][0]-x[i-2][0]+0.5*h*h*f(t[i-1],x[i-1])[1]
v_i=(x_i-x[i-1][0])/h
x[i]=np.array([x_i,v_i])
i+=1
return t,x
#tłumienie
#t,x=EE(lambda x:-x,1.0,0.0,100.0,0.01)
#py.plot(t,x)
#py.show()
#oscylator harmoniczny
def oscylator(t,y):
x=y[0]
xdot=y[1]
return np.array([xdot,-x])
#t,x=EE(oscylator,[1.0,1.0],0.0,1000.0,0.01)
t,x=leapfrog(oscylator,[1.0,1.0],0.0,100000.0,0.1)
py.plot(x[:,0],x[:,1])
py.show()
#Ec=map(lambda z:z[0]**2+z[1]**2,x)
#py.plot(t,Ec)
py.show()
Zadanie
Rozwiąż układ równań różniczkowych Lorenza dany wzorami
[math] \begin{cases}\dot x=\sigma y-\sigma x\\\dot y=-xz+rx-y\\\dot z=xy-bz\end{cases}, [/math]
metodą całkowania Rungego–Kutty drugiego rzędu z α = 2/3, czyli
[math] k_1 &= f(t_n,y_n), \\k_2 &= f(t_n + \tfrac{2}{3}h, y_n + \tfrac{2}{3}h k_1), \\y_{n+1} &= y_n + h\left(\tfrac{1}{4}k_1+\tfrac{3}{4}k_2\right). [/math]
Przyjmij sigma=10, b=8/3, r=99.96, krok czasowy h=0.005 i warunki początkowe x=1,y=0,z=0. Wykonaj 8000 kroków czasowych. Układ po pewnym czasie zacznie poruszać się po pewnej periodycznej trajektorii. Wykonaj 3 rysunki TEJ PERIODYCZNEJ TRAJEKTORII (bez okresu dochodzenia do niej) w płaszczyznach (x,y), (y,z) i (z,x). Wypisz na ekranie przedziały wartości jakie przyjmują zmienne x,y i z na periodycznej trajektorii oraz okres trajektorii periodycznej.
import numpy as np
import pylab as py
#rozwiazujemy rownanie dx(t)/dt=f(t,x)
#metoda Explicit Euler
#f - funkcja f z rownania
#x0-wartosc poczatkowa
#t0-czas poczatkowy
#tk-czas koncowy
#h-krok czasowy
def RK2(f,x0,t0,tk,h):
#generujemy wektor czasow
t=np.arange(t0,tk,h)
#liczba krokow czasowych
N=len(t)
#wektor wynikowy
if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
else: x=np.zeros(N)
#wpisujemy wartosc poczatkowa
x[0]=np.array(x0)
#index
i=1
#petla glowna
while (i<N):
k1=h*f(t[i-1],x[i-1])
k2=h*f(t[i-1]+h*2.0/3.0,x[i-1]+k1*2.0/3.0)
x[i]=np.array(x[i-1]+0.25*k1+0.75*k2)
i+=1
return t,x
#uklad Lorenza
def Lorenza(t,y):
sigma=10.0
b=8.0/3
r=99.96
xx=y[0]
yy=y[1]
zz=y[2]
xdot=sigma*(yy-xx)
ydot=-xx*zz+r*xx-yy
zdot=xx*yy-b*zz
return np.array([xdot,ydot,zdot])
t,x=RK2(Lorenza,[1.0,0.0,0.0],0.0,40.0,0.005)
py.plot(x[2500:,0],x[2500:,1])
py.show()
py.plot(x[2500:,1],x[2500:,2])
py.show()
py.plot(x[2500:,2],x[2500:,0])
py.show()
#szukanie okresu
prog=140
lista=[]
for i in range(2500,8000):
if (x[i-1,2]<prog) and (x[i,2]>prog): lista.append(i)
print 'okres to:',np.mean(np.diff(lista)*0.005)
print 'zmienna x przyjmuje wartosci z zakresu: (',min(x[2500:,0]),',',max(x[2500:,0]),')'
print 'zmienna y przyjmuje wartosci z zakresu: (',min(x[2500:,1]),',',max(x[2500:,1]),')'
print 'zmienna z przyjmuje wartosci z zakresu: (',min(x[2500:,2]),',',max(x[2500:,2]),')'
