TI/Programowanie dla Fizyków Medycznych/RRZ: Różnice pomiędzy wersjami

Z Brain-wiki
Linia 145: Linia 145:
 
py.show()
 
py.show()
 
</source>
 
</source>
 +
===Wahadło===
 +
<source lang="python">
 +
import numpy as np
 +
import pylab as py
 +
 +
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
 +
 +
 +
#oscylator harmoniczny
 +
def oscylator(t,y):
 +
    f0=1.0
 +
    w0=1.0
 +
   
 +
    x=y[0]
 +
    xdot=y[1]
 +
    return np.array([xdot,f0*np.cos(oscylator.W*t)-oscylator.Gamma*xdot-w0*w0*x])
 +
 +
def amplituda(x):
 +
    lista=x[2000:,0]
 +
    return (max(lista)-min(lista))*0.5
 +
 +
def okres(t,x):
 +
    t=t[2000:]
 +
    x=x[2000:,0]
 +
    czasyMaksimow=[]
 +
    for i in xrange(1,len(x)-1):
 +
        if x[i]>x[i-1] and x[i]>x[i+1]:
 +
            czasyMaksimow.append(t[i])
 +
    return np.mean(np.diff(czasyMaksimow))
 +
   
 +
 +
oscylator.W=1.0
 +
oscylator.Gamma=0.1
 +
 +
#Amplituda wahań w funkcji Omega (częstości wymuszania)
 +
#Omegas=np.arange(0.1,3.0,0.05)
 +
#amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]
 +
#py.plot(Omegas,amp)
 +
#py.show()
 +
 +
#Okres wahań w funkcji Omega (częstości wymuszania)
 +
#Omegas=np.arange(0.1,3.0,0.1)
 +
#okresy=[okres(*RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)) for oscylator.W in Omegas]
 +
#py.plot(Omegas,okresy)
 +
#py.show()
 +
 +
#Amplituda wahań w funkcji wsp. tłumienia Gamma
 +
Gammas=np.arange(0.1,3.0,0.1)
 +
amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.Gamma in Gammas]
 +
py.plot(Gammas,amp)
 +
py.show()
 +
</source>
 +
 +
  
 
===Zadanie===
 
===Zadanie===

Wersja z 21:13, 3 cze 2015

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()

Wahadło

import numpy as np
import pylab as py

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


#oscylator harmoniczny
def oscylator(t,y):
    f0=1.0
    w0=1.0
    
    x=y[0]
    xdot=y[1]
    return np.array([xdot,f0*np.cos(oscylator.W*t)-oscylator.Gamma*xdot-w0*w0*x])

def amplituda(x):
    lista=x[2000:,0]
    return (max(lista)-min(lista))*0.5

def okres(t,x):
    t=t[2000:]
    x=x[2000:,0]
    czasyMaksimow=[]
    for i in xrange(1,len(x)-1):
        if x[i]>x[i-1] and x[i]>x[i+1]:
            czasyMaksimow.append(t[i])
    return np.mean(np.diff(czasyMaksimow))
    

oscylator.W=1.0
oscylator.Gamma=0.1

#Amplituda wahań w funkcji Omega (częstości wymuszania)
#Omegas=np.arange(0.1,3.0,0.05)
#amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]
#py.plot(Omegas,amp)
#py.show()

#Okres wahań w funkcji Omega (częstości wymuszania)
#Omegas=np.arange(0.1,3.0,0.1)
#okresy=[okres(*RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)) for oscylator.W in Omegas]
#py.plot(Omegas,okresy)
#py.show()

#Amplituda wahań w funkcji wsp. tłumienia Gamma
Gammas=np.arange(0.1,3.0,0.1)
amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.Gamma in Gammas]
py.plot(Gammas,amp)
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]),')'