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

Z Brain-wiki
 
(Nie pokazano 32 wersji utworzonych przez 3 użytkowników)
Linia 1: Linia 1:
 +
__NOTOC__
 
==Równania różniczkowe zwyczajne==
 
==Równania różniczkowe zwyczajne==
Zajmieny się teraz problemem numerycznego rozwiązywania równań różniczkowych zwyczajnych o postaci:
+
Zajmiemy się teraz problemem numerycznego rozwiązywania równań różniczkowych zwyczajnych o postaci:
 +
 
 +
<!--
 +
{{Hide in print|This text will not be shown in the print.}}
 +
{{Only in print|This text will only be shown in the print.}}
 +
-->
  
 
<math> \frac{dy(t)}{dt}  = f(t,y(t))</math>,
 
<math> \frac{dy(t)}{dt}  = f(t,y(t))</math>,
Linia 16: Linia 22:
 
<math> \frac{d}{dt} \binom{x(t)}{x'(t)}  = \binom{x'(t)}{\omega(t,x(t))}</math>.
 
<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,
+
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 pochodną 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łcając uzyskujemy
  
<source lang="python">
+
<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}  = 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
 +
 
 +
<!--<source lang="python">-->
 +
<syntaxhighlight lang="python">
 
import numpy as np
 
import numpy as np
 
import pylab as py
 
import pylab as py
Linia 31: Linia 57:
 
#h-krok czasowy
 
#h-krok czasowy
  
def EE(f,x0,t0,tk,h):  
+
def EE(f,x0,t0,tk,h):
 
     #generujemy wektor czasow
 
     #generujemy wektor czasow
 
     t=np.arange(t0,tk,h)
 
     t=np.arange(t0,tk,h)
 +
 
     #liczba krokow czasowych
 
     #liczba krokow czasowych
 
     N=len(t)
 
     N=len(t)
 +
 
     #wektor wynikowy
 
     #wektor wynikowy
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
+
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0))) # gdy mamy do czynienia w równaniem wektorowym
     else: x=np.zeros(N)
+
     else: x=np.zeros(N) #dla przypadku skalarnego
 +
 
 
     #wpisujemy wartosc poczatkowa
 
     #wpisujemy wartosc poczatkowa
 
     x[0]=np.array(x0)
 
     x[0]=np.array(x0)
Linia 48: Linia 77:
 
         i+=1
 
         i+=1
 
     return t,x
 
     return t,x
 +
</syntaxhighlight>
 +
<!--</source>-->
 +
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
 +
<source lang="python">
 +
def oscylator(t,y):
 +
    x=y[0]
 +
    xdot=y[1]
 +
    return np.array([xdot,-x])
 +
</source>
 +
Rozwiążmy to równanie z warunkiem 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 następująco
 +
 +
[[Plik:img01.png]]
 +
 +
Jeżeli zaś wydłużymy czas symulacji do 1000 otrzymamy
 +
 +
[[Plik:img2.png]]
 +
 +
Amplituda oscylacji rośnie wykładniczo i rozwiązanie numeryczne bardzo szybko przestaje mieć cokolwiek wspólnego ze ścisłym rozwiązaniem, którego amplituda jest przecież stała. Metoda Explicit Euler już po kilku krokach czasowych przestaje przypominać ścisłe rozwiązanie. Niestety trudno jest zupełnie wyeliminować to zjawisko, za to możemy użyć metody, która znacznie wolniej będzie się oddalać od ścisłego rozwiązania. Zauważmy, że w metodzie Explicit Euler w każdym kroku czasowym tylko raz liczyliśmy wartość funkcji f. Liczbę wywołań funkcji f w każdym kroku czasowym nazywamy rzędem metody, stąd Explicit Euler jest metodą pierwszego rzędu. Wprowadźmy teraz przykładowe metody rzędu drugiego.
 +
 +
==Metoda Żabiego Skoku== 
 +
W poprzedniej metodzie liczyliśmy wartość funkcji f w chwili <math>t_n</math>, która była pochodną po czasie naszego ścisłego rozwiązania. Kolejny punkt <math>x_{n+1}</math> był liczony z przybliżenia liniowego funkcji  w chwili poprzedniej. Jeżeli faktyczna trajektoria ma niezerową drugą pochodną to takie liniowe przybliżenie zawsze będzie nas oddalało od ścisłego rozwiązania. Dosyć prostym pomysłem na poprawienie zbieżności metody jest tak zwany żabi skok. Policzmy najpierw wartość zmiennej x przesuwając się w czasie o h/2 i policzmy wówczas pochodną, którą oznaczmy przez <math>k_2</math>
 +
 +
<math> k_1=f(t_n,x_n) </math>.
 +
 +
<math> k_2=f(t_n+h/2,x_n+h/2*k_1) </math>.
  
 +
Następnie używamy pochodnej <math> k_2 </math> zamiast pochodnej <math> k_1 </math> do obliczenia wartości funkcji w kolejnym kroku czasowym.
 +
 +
<math> x_{n+1}  = x_n + h k_2 </math>.
 +
 +
Przykładowa implementacja tej metody wygląda następująco
 +
<source lang="python">
 
def leapfrog(f,x0,t0,tk,h):  
 
def leapfrog(f,x0,t0,tk,h):  
 +
 
     #generujemy wektor czasow
 
     #generujemy wektor czasow
 
     t=np.arange(t0,tk,h)
 
     t=np.arange(t0,tk,h)
 +
 
     #liczba krokow czasowych
 
     #liczba krokow czasowych
 
     N=len(t)
 
     N=len(t)
 +
 
     #wektor wynikowy
 
     #wektor wynikowy
 
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
 
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
 
     else: x=np.zeros(N)
 
     else: x=np.zeros(N)
 +
 
     #wpisujemy wartosc poczatkowa
 
     #wpisujemy wartosc poczatkowa
 
     x[0]=np.array(x0)
 
     x[0]=np.array(x0)
 +
 
     #index
 
     #index
 
     i=1
 
     i=1
 +
 
     #petla glowna
 
     #petla glowna
 
     while (i<N):
 
     while (i<N):
Linia 68: Linia 140:
 
         i+=1
 
         i+=1
 
     return t,x
 
     return t,x
 +
</source>
 +
Rozwiązanie równania oscylatora tą metodą dla identycznych jak poprzednio czasów da następujące wyniki
 +
<source lang="python">
 +
t,x=leapfrog(oscylator,[1.0,1.0],0.0,100,0.01)
 +
py.plot(t,x[:,0])
 +
py.show()
 +
</source>
 +
 +
[[Plik:img3.png]]
 +
 +
<source lang="python">
 +
t,x=leapfrog(oscylator,[1.0,1.0],0.0,1000,0.01)
 +
py.plot(t,x[:,0])
 +
py.show()
 +
</source>
 +
 +
[[Plik:img4.png]]
 +
 +
==Metoda Heuna==
 +
Kolejną metodą niewiele różniącą się od poprzedniej jest metoda Heuna. Zdefiniowana jest ona przez równania
 +
<math> k_1=f(t_n,x_n) </math>,
 +
 +
<math> k_2=f(t_n+h/2,x_n+h/2*k_1) </math>,
 +
 +
<math> x_{n+1}  = x_n + h k_2 </math>.
 +
 +
Implementacja wygląda następująco
 +
 +
<source lang="python">
  
 
def Heun(f,x0,t0,tk,h):  
 
def Heun(f,x0,t0,tk,h):  
Linia 88: Linia 189:
 
         i+=1
 
         i+=1
 
     return t,x
 
     return t,x
 +
</source>
  
def RK4(f,x0,t0,tk,h):
+
==Runge-Kutta czwartego rzędu==
    #generujemy wektor czasow
+
Ostatnią metodą, którą omówimy jest najbardziej popularna metoda zwana w skrócie RK4. Metoda ta uznawana jest za kanoniczną i w większości zastosowań dającą najlepsze wyniki. Metody wyższego rzędu nie wnoszą już do wyniku znaczącej poprawy. Jak sugeruje nazwa metody, jej rząd to 4, czyli w każdym kroku czasowym czterokrotnie wywołujemy funkcję f. Metoda ta zdefiniowana jest wzorami
    t=np.arange(t0,tk,h)
+
 
    #liczba krokow czasowych
+
<math>    k_1 = f \left( t_n, x_n \right) </math>,
    N=len(t)
+
 
    #wektor wynikowy
+
<math>    k_2 = f \left( t_n + {h \over 2}, x_n + {1 \over 2} k_1 \right) </math>,
    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:
+
<math>    k_3 = f \left( t_n + {h \over 2}, x_n + {1 \over 2} k_2 \right) </math>,
# 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
+
<math>    k_4 = f \left( t_n + h, x_n + k_3 \right) </math>,
#t,x=EE(lambda x:-x,1.0,0.0,100.0,0.01)
 
#py.plot(t,x)
 
#py.show()
 
  
 +
<math>    x_{n+1} = x_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4) </math>,
  
#oscylator harmoniczny
+
a implementacja wygląda następująco
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()
 
</source>
 
===Wahadło===
 
 
<source lang="python">
 
<source lang="python">
import numpy as np
 
import pylab as py
 
 
 
def RK4(f,x0,t0,tk,h):  
 
def RK4(f,x0,t0,tk,h):  
 
     #generujemy wektor czasow
 
     #generujemy wektor czasow
Linia 190: Linia 228:
 
         i+=1
 
         i+=1
 
     return t,x
 
     return t,x
 +
</source>
  
 +
==Przykłady==
 +
===Zadanie - Wahadło matematyczne z tłumieniem i siłą wymuszającą===
  
#oscylator harmoniczny
+
Rozwiąż numerycznie metodą RK4 równanie różniczkowe oscylatora harmonicznego z tłumieniem i siłą wymuszającą
 +
 
 +
<math> \frac{d^2x}{dt^2} + \Gamma \frac{dx}{dt} + w_0^2 x = f_0 \cos(W t) </math>,
 +
 
 +
przyjmując parametry <math> f_0 =1, w_0=1, \Gamma=0.1, h=0.1 </math>. Wykreśl zależność amplitudy drgań w funkcji częstości siły wymuszającej W, dla W z przedziału [0.1,3].
 +
===Rozwiązanie===
 +
Zacznijmy od sprowadzenia równania drugiego stopnia do równania pierwszego stopnia i zapisania go w postaci funkcji
 +
 
 +
<source lang="python">
 
def oscylator(t,y):
 
def oscylator(t,y):
 
     f0=1.0
 
     f0=1.0
 
     w0=1.0
 
     w0=1.0
      
+
     Gamma=0.1
 +
 
 
     x=y[0]
 
     x=y[0]
 
     xdot=y[1]
 
     xdot=y[1]
 
     return np.array([xdot,f0*np.cos(oscylator.W*t)-oscylator.Gamma*xdot-w0*w0*x])
 
     return np.array([xdot,f0*np.cos(oscylator.W*t)-oscylator.Gamma*xdot-w0*w0*x])
  
 +
oscylator.W=1.0
 +
</source>
 +
Nieprzypadkowo parametr W nie jest definiowany w samej funkcji jako zmienna wewnętrzna, ale jako atrybut obiektu jakim jest funkcja. Dzięki takiej konstrukcji łatwo będzie nam zmieniać parametr W, czyli częstość siły wymuszającej. Zobaczmy jak wygląda trajektoria będąca rozwiązaniem tego równania dla warunku początkowego [0,1] i czasu końcowego równego 400.
 +
<source lang="python">
 +
py.plot(*RK4(oscylator,[0.0,1.0],0.0,400.0,0.1))
 +
py.show()
 +
</source>
 +
 +
[[Plik:oscy1.png]]
 +
 +
Widać, że początkowo układ dochodzi do stanu regularnych oscylacji. W analizie amplitudy interesuje nas jedynie końcowa część więc ograniczymy się do analizy trajektorii od czasu 200 do czasu 400. Napiszmy teraz funkcję, która na podstawie trajektorii wyznaczy nam amplitudę oscylacji
 +
 +
<source lang="python">
 
def amplituda(x):
 
def amplituda(x):
 
     lista=x[2000:,0]
 
     lista=x[2000:,0]
 
     return (max(lista)-min(lista))*0.5
 
     return (max(lista)-min(lista))*0.5
 +
</source>
  
def okres(t,x):
+
Interesować nas będzie amplituda drgań w funkcji częstość W. Wygenerujmy listę wartości W, dla których będziemy liczyć amplitudę.
    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
+
<source lang="python">
oscylator.Gamma=0.1
+
Omegas=np.arange(0.1,3.0,0.05)
 
+
</source>
#Amplituda wahań w funkcji Omega (częstości wymuszania)
+
Możemy teraz dla każdej z wartości W rozwiązać numerycznie równanie różniczkowe i wyznaczyć odpowiadającą amplitudę oscylacji
#Omegas=np.arange(0.1,3.0,0.05)
+
<source lang="python">
#amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]
+
amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]
#py.plot(Omegas,amp)
+
</source>
#py.show()
+
Wynik końcowy wygląda następująco:
 
+
<source lang="python">
#Okres wahań w funkcji Omega (częstości wymuszania)
+
py.plot(Omegas,amp)
#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()
 
py.show()
 
</source>
 
</source>
 +
[[Plik:oscy2.png]]
  
 +
Jak można było się domyślić amplituda jest największa gdy częstotliwość wymuszania W pokrywa się z wartością częstotliwości drgań własnych <math> w_0=1 </math>
  
 
+
===Zadanie - Układ Lorenza===
===Zadanie===
 
 
Rozwiąż układ równań różniczkowych Lorenza dany wzorami
 
Rozwiąż układ równań różniczkowych Lorenza dany wzorami
  
Linia 246: Linia 294:
 
metodą całkowania Rungego–Kutty drugiego rzędu z α = 2/3, czyli
 
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>  
+
<math>k_1 = f(t_n,x_n)</math> ,  
 +
 
 +
<math>k_2 = f(t_n + \tfrac{2}{3}h, x_n + \tfrac{2}{3}h k_1)</math> ,  
 +
 
 +
<math>x_{n+1} = x_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.
 
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.
 
+
===Rozwiązanie===
 +
Zacznijmy od implementacji podanej w treści zadania metody całkowania RK2
 
<source lang="python">
 
<source lang="python">
 
import numpy as np
 
import numpy as np
Linia 264: Linia 317:
  
 
def RK2(f,x0,t0,tk,h):  
 
def RK2(f,x0,t0,tk,h):  
 +
 
     #generujemy wektor czasow
 
     #generujemy wektor czasow
 
     t=np.arange(t0,tk,h)
 
     t=np.arange(t0,tk,h)
 +
 
     #liczba krokow czasowych
 
     #liczba krokow czasowych
 
     N=len(t)
 
     N=len(t)
 +
 
     #wektor wynikowy
 
     #wektor wynikowy
 
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
 
     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0)))
 
     else: x=np.zeros(N)
 
     else: x=np.zeros(N)
 +
 
     #wpisujemy wartosc poczatkowa
 
     #wpisujemy wartosc poczatkowa
 
     x[0]=np.array(x0)
 
     x[0]=np.array(x0)
 +
 
     #index
 
     #index
 
     i=1
 
     i=1
 +
 
     #petla glowna
 
     #petla glowna
 
     while (i<N):
 
     while (i<N):
Linia 282: Linia 341:
 
         i+=1
 
         i+=1
 
     return t,x
 
     return t,x
 
+
</source>
#uklad Lorenza
+
Następnie napiszmy funkcję opisującą układ Lorenza
 +
<source lang="python">
 
def Lorenza(t,y):
 
def Lorenza(t,y):
 
     sigma=10.0
 
     sigma=10.0
Linia 295: Linia 355:
 
     zdot=xx*yy-b*zz
 
     zdot=xx*yy-b*zz
 
     return np.array([xdot,ydot,zdot])
 
     return np.array([xdot,ydot,zdot])
 +
</source>
  
 +
Zobaczmy teraz jak wyglądają trajektorie wszystkich trzech współrzędnych  rozwiązania z zadanymi parametrami
 +
<source lang="python">
 
t,x=RK2(Lorenza,[1.0,0.0,0.0],0.0,40.0,0.005)
 
t,x=RK2(Lorenza,[1.0,0.0,0.0],0.0,40.0,0.005)
 +
py.plot(t,x[:,1])
 +
py.show()
 +
py.plot(t,x[:,2])
 +
py.show()
 +
py.plot(t,x[:,0])
 +
py.show()
 +
</source>
 +
 +
[[Plik:lor1.png]]
 +
 +
[[Plik:lor2.png]]
 +
 +
[[Plik:lor3.png]]
 +
 +
Możemy zauważyć, że układ początkowo zachowuje się chaotycznie, a potem dąży do pewnego stanu ustalonego (tzw. atraktora). Cała trajektoria składa się z 8000 punktów, przyjmijmy że powyżej punktu o numerze 2500 mamy już do czynienia tylko z periodyczną trajektorią. Wykreślmy zatem portrety fazowe, o których mowa w treści zadania
 +
<source lang="python">
 
py.plot(x[2500:,0],x[2500:,1])
 
py.plot(x[2500:,0],x[2500:,1])
 
py.show()
 
py.show()
Linia 303: Linia 382:
 
py.plot(x[2500:,2],x[2500:,0])
 
py.plot(x[2500:,2],x[2500:,0])
 
py.show()
 
py.show()
 +
</source>
 +
 +
[[Plik:lor4.png]]
 +
 +
[[Plik:lor5.png]]
 +
 +
[[Plik:lor6.png]]
  
#szukanie okresu
+
Wypisanie zakresów jest już tylko formalnością
 +
<source lang="python">
 +
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]),')'
 +
</source>
 +
Okres trajektorii periodycznej możemy znaleźć na przykład w ten sposób
 +
<source lang="python">
 
prog=140
 
prog=140
 
lista=[]
 
lista=[]
Linia 310: Linia 403:
 
     if (x[i-1,2]<prog) and (x[i,2]>prog): lista.append(i)
 
     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 '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]),')'
 
  
 +
>>> okres to: 1.0975
 +
>>> zmienna x przyjmuje wartosci z zakresu: ( -33.4431203059 , 25.8037953495 )
 +
>>> zmienna y przyjmuje wartosci z zakresu: ( -56.7169238157 , 37.3166709986 )
 +
>>> zmienna z przyjmuje wartosci z zakresu: ( 53.4652816712 , 144.264397579 )
  
 
</source>
 
</source>
  
 
[["Programowanie dla Fizyków Medycznych"]]
 
[["Programowanie dla Fizyków Medycznych"]]

Aktualna wersja na dzień 14:04, 12 cze 2015

Równania różniczkowe zwyczajne

Zajmiemy 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 pochodną 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łcają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} = 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 do czynienia 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 warunkiem 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 następująco

Img01.png

Jeżeli zaś wydłużymy czas symulacji do 1000 otrzymamy

Img2.png

Amplituda oscylacji rośnie wykładniczo i rozwiązanie numeryczne bardzo szybko przestaje mieć cokolwiek wspólnego ze ścisłym rozwiązaniem, którego amplituda jest przecież stała. Metoda Explicit Euler już po kilku krokach czasowych przestaje przypominać ścisłe rozwiązanie. Niestety trudno jest zupełnie wyeliminować to zjawisko, za to możemy użyć metody, która znacznie wolniej będzie się oddalać od ścisłego rozwiązania. Zauważmy, że w metodzie Explicit Euler w każdym kroku czasowym tylko raz liczyliśmy wartość funkcji f. Liczbę wywołań funkcji f w każdym kroku czasowym nazywamy rzędem metody, stąd Explicit Euler jest metodą pierwszego rzędu. Wprowadźmy teraz przykładowe metody rzędu drugiego.

Metoda Żabiego Skoku

W poprzedniej metodzie liczyliśmy wartość funkcji f w chwili [math]t_n[/math], która była pochodną po czasie naszego ścisłego rozwiązania. Kolejny punkt [math]x_{n+1}[/math] był liczony z przybliżenia liniowego funkcji w chwili poprzedniej. Jeżeli faktyczna trajektoria ma niezerową drugą pochodną to takie liniowe przybliżenie zawsze będzie nas oddalało od ścisłego rozwiązania. Dosyć prostym pomysłem na poprawienie zbieżności metody jest tak zwany żabi skok. Policzmy najpierw wartość zmiennej x przesuwając się w czasie o h/2 i policzmy wówczas pochodną, którą oznaczmy przez [math]k_2[/math]

[math] k_1=f(t_n,x_n) [/math].

[math] k_2=f(t_n+h/2,x_n+h/2*k_1) [/math].

Następnie używamy pochodnej [math] k_2 [/math] zamiast pochodnej [math] k_1 [/math] do obliczenia wartości funkcji w kolejnym kroku czasowym.

[math] x_{n+1} = x_n + h k_2 [/math].

Przykładowa implementacja tej metody wygląda następują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

Rozwiązanie równania oscylatora tą metodą dla identycznych jak poprzednio czasów da następujące wyniki

t,x=leapfrog(oscylator,[1.0,1.0],0.0,100,0.01)
py.plot(t,x[:,0])
py.show()

Img3.png

t,x=leapfrog(oscylator,[1.0,1.0],0.0,1000,0.01)
py.plot(t,x[:,0])
py.show()

Img4.png

Metoda Heuna

Kolejną metodą niewiele różniącą się od poprzedniej jest metoda Heuna. Zdefiniowana jest ona przez równania [math] k_1=f(t_n,x_n) [/math],

[math] k_2=f(t_n+h/2,x_n+h/2*k_1) [/math],

[math] x_{n+1} = x_n + h k_2 [/math].

Implementacja wygląda następująco

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

Runge-Kutta czwartego rzędu

Ostatnią metodą, którą omówimy jest najbardziej popularna metoda zwana w skrócie RK4. Metoda ta uznawana jest za kanoniczną i w większości zastosowań dającą najlepsze wyniki. Metody wyższego rzędu nie wnoszą już do wyniku znaczącej poprawy. Jak sugeruje nazwa metody, jej rząd to 4, czyli w każdym kroku czasowym czterokrotnie wywołujemy funkcję f. Metoda ta zdefiniowana jest wzorami

[math] k_1 = f \left( t_n, x_n \right) [/math],

[math] k_2 = f \left( t_n + {h \over 2}, x_n + {1 \over 2} k_1 \right) [/math],

[math] k_3 = f \left( t_n + {h \over 2}, x_n + {1 \over 2} k_2 \right) [/math],

[math] k_4 = f \left( t_n + h, x_n + k_3 \right) [/math],

[math] x_{n+1} = x_n + {h \over 6} (k_1 + 2k_2 + 2k_3 + k_4) [/math],

a implementacja wygląda następująco

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

Przykłady

Zadanie - Wahadło matematyczne z tłumieniem i siłą wymuszającą

Rozwiąż numerycznie metodą RK4 równanie różniczkowe oscylatora harmonicznego z tłumieniem i siłą wymuszającą

[math] \frac{d^2x}{dt^2} + \Gamma \frac{dx}{dt} + w_0^2 x = f_0 \cos(W t) [/math],

przyjmując parametry [math] f_0 =1, w_0=1, \Gamma=0.1, h=0.1 [/math]. Wykreśl zależność amplitudy drgań w funkcji częstości siły wymuszającej W, dla W z przedziału [0.1,3].

Rozwiązanie

Zacznijmy od sprowadzenia równania drugiego stopnia do równania pierwszego stopnia i zapisania go w postaci funkcji

def oscylator(t,y):
    f0=1.0
    w0=1.0
    Gamma=0.1

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

oscylator.W=1.0

Nieprzypadkowo parametr W nie jest definiowany w samej funkcji jako zmienna wewnętrzna, ale jako atrybut obiektu jakim jest funkcja. Dzięki takiej konstrukcji łatwo będzie nam zmieniać parametr W, czyli częstość siły wymuszającej. Zobaczmy jak wygląda trajektoria będąca rozwiązaniem tego równania dla warunku początkowego [0,1] i czasu końcowego równego 400.

py.plot(*RK4(oscylator,[0.0,1.0],0.0,400.0,0.1))
py.show()

Oscy1.png

Widać, że początkowo układ dochodzi do stanu regularnych oscylacji. W analizie amplitudy interesuje nas jedynie końcowa część więc ograniczymy się do analizy trajektorii od czasu 200 do czasu 400. Napiszmy teraz funkcję, która na podstawie trajektorii wyznaczy nam amplitudę oscylacji

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

Interesować nas będzie amplituda drgań w funkcji częstość W. Wygenerujmy listę wartości W, dla których będziemy liczyć amplitudę.

Omegas=np.arange(0.1,3.0,0.05)

Możemy teraz dla każdej z wartości W rozwiązać numerycznie równanie różniczkowe i wyznaczyć odpowiadającą amplitudę oscylacji

amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]

Wynik końcowy wygląda następująco:

py.plot(Omegas,amp)
py.show()

Oscy2.png

Jak można było się domyślić amplituda jest największa gdy częstotliwość wymuszania W pokrywa się z wartością częstotliwości drgań własnych [math] w_0=1 [/math]

Zadanie - Układ Lorenza

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,x_n)[/math] ,

[math]k_2 = f(t_n + \tfrac{2}{3}h, x_n + \tfrac{2}{3}h k_1)[/math] ,

[math]x_{n+1} = x_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.

Rozwiązanie

Zacznijmy od implementacji podanej w treści zadania metody całkowania RK2

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

Następnie napiszmy funkcję opisującą układ 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])

Zobaczmy teraz jak wyglądają trajektorie wszystkich trzech współrzędnych rozwiązania z zadanymi parametrami

t,x=RK2(Lorenza,[1.0,0.0,0.0],0.0,40.0,0.005)
py.plot(t,x[:,1])
py.show()
py.plot(t,x[:,2])
py.show()
py.plot(t,x[:,0])
py.show()

Lor1.png

Lor2.png

Lor3.png

Możemy zauważyć, że układ początkowo zachowuje się chaotycznie, a potem dąży do pewnego stanu ustalonego (tzw. atraktora). Cała trajektoria składa się z 8000 punktów, przyjmijmy że powyżej punktu o numerze 2500 mamy już do czynienia tylko z periodyczną trajektorią. Wykreślmy zatem portrety fazowe, o których mowa w treści zadania

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

Lor4.png

Lor5.png

Lor6.png

Wypisanie zakresów jest już tylko formalnością

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]),')'

Okres trajektorii periodycznej możemy znaleźć na przykład w ten sposób

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)

>>> okres to: 1.0975
>>> zmienna x przyjmuje wartosci z zakresu: ( -33.4431203059 , 25.8037953495 )
>>> zmienna y przyjmuje wartosci z zakresu: ( -56.7169238157 , 37.3166709986 )
>>> zmienna z przyjmuje wartosci z zakresu: ( 53.4652816712 , 144.264397579 )

"Programowanie dla Fizyków Medycznych"