TI/Programowanie dla Fizyków Medycznych/RRZ: Różnice pomiędzy wersjami
Linia 78: | Linia 78: | ||
return np.array([xdot,-x]) | return np.array([xdot,-x]) | ||
</source> | </source> | ||
+ | Rozwiążmy to równanie z warunkeim początkowym [1.0,1.0] i od czasu od 0 do 100. | ||
+ | <source lang="python"> | ||
+ | t,x=EE(oscylator,[1.0,1.0],0.0,100,0.01) | ||
+ | py.plot(t,x[:,0]) | ||
+ | py.show() | ||
+ | </source> | ||
+ | rozwiązanie wygląda wówczas nastepująco | ||
+ | [[Plik:img01.png]] | ||
− | |||
<source lang="python"> | <source lang="python"> | ||
def leapfrog(f,x0,t0,tk,h): | def leapfrog(f,x0,t0,tk,h): |
Wersja z 13:34, 5 cze 2015
Równania różniczkowe zwyczajne
Zajmieny się teraz problemem numerycznego rozwiązywania równań różniczkowych zwyczajnych o postaci:
[math] \frac{dy(t)}{dt} = f(t,y(t))[/math],
z warunkeim początkowym
[math]y(t_0)=y_0[/math].
Zauważmy, że przykładowe równanie różniczkowe drugiego rzędu
[math] \frac{d^2 x(t)}{dt^2} = \omega(t,x(t))[/math],
można zapisać jako
[math] \frac{d}{dt} \binom{x(t)}{x'(t)} = \binom{x'(t)}{\omega(t,x(t))}[/math].
W analogiczny sposób równanie dowolnego rzędy możemy zapisać jako wektorowe równanie różniczkowe pierwszego rzędu. Wystarczy zatem, że skupimy się na rozwiązywaniu równań pierwszego rzędu, Rozwiązaniem postawionego problemu są ciągłe funkcje zmiennej czasowej t. Rozwiązanie numeryczne takiego problemu ogranicza się jednak do znalezienia wartości funkcji y(t) w skończonej liczbie punktów czasowych. W najprostrzym przypadku (do którego się tutaj ograniczymy) zakładamy, że punkty te są od siebie równo oddalone, a odległość między nimi nazywamy krokiem czasowym i tradycyjnie oznaczamy literą h. Zatem rozwiązanie równania na przedziale [math](t_0,t_k)[/math] sprwadzamy do rozwiązania w sekwencji czasów [math]t_0, t+1=t_0+h,t_2=t_0+2h,...,t_k=Nh[/math]. Poprzez [math]x_n[/math] oznaczać będziemy numeryczne przybliżenie ścisłego rozwiązania [math]x(t_n)[/math].
Metoda Eulera
Najprostszą metodą numeryczną rozwiązywania równań różniczkowych jest metoda Eulera. Przybliżmy pochodzną czasową występującą po lewej stronie równania przez iloraz różnicowy
[math] \frac{d x(t)}{dt} \approx \frac{x(t+h)-x(t)}{h}[/math],
przekształacjąc uzyskujemy
[math] x(t+h) \approx x(t)+ h \frac{d x(t)}{dt} [/math]
a po podstawieniu rozwiązywanego równania mamy
[math] x(t+h) \approx x(t)+ h f(t,x(t)) [/math].
Możemy to zapisać w postaci dyskretnej
[math] x_{n+1} \approx x_n + h f(t_n,x_n) [/math].
Wartość w kolejnej chwili czasu dana jest explicite poprzez wartość w chwili poprzedniej. Metoda ta nazywa się Explicit Euler. Możemy teraz zaimplementować ją w pythonie
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))) # gdy mamy doczynienia w równaniem wektorowym
else: x=np.zeros(N) #dla przypadku skalarnego
#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
Najłatwiej będzie przetestować napisaną metodę na równaniu którego ścisłe rozwiązanie jest znane.Zacznijmy zatem od równania oscylatora harmonicznego
def oscylator(t,y):
x=y[0]
xdot=y[1]
return np.array([xdot,-x])
Rozwiążmy to równanie z warunkeim początkowym [1.0,1.0] i od czasu od 0 do 100.
t,x=EE(oscylator,[1.0,1.0],0.0,100,0.01)
py.plot(t,x[:,0])
py.show()
rozwiązanie wygląda wówczas nastepująco
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()
#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]),')'