TI/Programowanie dla Fizyków Medycznych/RRZ - Historia wersji

SuperAdmin: /* Metoda Eulera */

2015-06-12T14:04:35Z

Metoda Eulera

SuperAdmin: /* Metoda Eulera */

2015-06-12T14:03:23Z

Metoda Eulera

Tgubiec: /* Rozwiązanie */

2015-06-09T12:59:35Z

Rozwiązanie

Tgubiec: /* Rozwiązanie */

2015-06-09T12:12:11Z

Rozwiązanie

Tgubiec: /* Rozwiązanie */

2015-06-09T12:09:21Z

Rozwiązanie

Tgubiec: /* Rozwiązanie */

2015-06-09T12:06:34Z

Rozwiązanie

Tgubiec: /* Runge-Kutta czwartego rzędu */

2015-06-09T12:03:57Z

Runge-Kutta czwartego rzędu

Tgubiec: /* Metoda Eulera */

2015-06-09T12:02:36Z

Metoda Eulera

Tgubiec: /* Metoda Żabiego Skoku */

2015-06-09T12:02:10Z

Metoda Żabiego Skoku

Tgubiec: /* Metoda Eulera */

2015-06-09T11:58:53Z

Metoda Eulera

@@ Linia 45: / Linia 45: @@
 <!--<source lang="python">-->
 <syntaxhighlight lang="python">
-import numpy as np<br>
+import numpy as np
-import pylab as py<br>
+import pylab as py
-<br>
-#rozwiazujemy rownanie dx(t)/dt=f(t,x)<br>
+#rozwiazujemy rownanie dx(t)/dt=f(t,x)
-<br>
-#metoda Explicit Euler<br>
+#metoda Explicit Euler
-#f - funkcja f z rownania<br>
+#f - funkcja f z rownania
-#x0-wartosc poczatkowa<br>
+#x0-wartosc poczatkowa
-#t0-czas poczatkowy<br>
+#t0-czas poczatkowy
-#tk-czas koncowy<br>
+#tk-czas koncowy
-#h-krok czasowy<br>
+#h-krok czasowy
-<br>
-def EE(f,x0,t0,tk,h):<br>
+def EE(f,x0,t0,tk,h):
-     #generujemy wektor czasow<br>
+     #generujemy wektor czasow
-     t=np.arange(t0,tk,h)<br>
+     t=np.arange(t0,tk,h)
-<br>
-     #liczba krokow czasowych<br>
+     #liczba krokow czasowych
-     N=len(t)<br>
+     N=len(t)
-<br>
-     #wektor wynikowy<br>
+     #wektor wynikowy
-     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0))) # gdy mamy do czynienia w równaniem wektorowym<br>
+     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<br>
+     else: x=np.zeros(N) #dla przypadku skalarnego
-<br>
-     #wpisujemy wartosc poczatkowa<br>
+     #wpisujemy wartosc poczatkowa
-     x[0]=np.array(x0)<br>
+     x[0]=np.array(x0)
-     #index<br>
+     #index
-     i=1<br>
+     i=1
-     #petla glowna<br>
+     #petla glowna
-     while (i<N):<br>
+     while (i<N):
-         x[i]=np.array(x[i-1]+h*f(t[i-1],x[i-1]))<br>
+         x[i]=np.array(x[i-1]+h*f(t[i-1],x[i-1]))
-         i+=1<br>
+         i+=1
-     return t,x<br>
+     return t,x
 </syntaxhighlight>
 <!--</source>-->

@@ Linia 45: / Linia 45: @@
 <!--<source lang="python">-->
 <syntaxhighlight lang="python">
-import numpy as np
+import numpy as np<br>
-import pylab as py
+import pylab as py<br>
+<br>
-#rozwiazujemy rownanie dx(t)/dt=f(t,x)
+#rozwiazujemy rownanie dx(t)/dt=f(t,x)<br>
+<br>
-#metoda Explicit Euler
+#metoda Explicit Euler<br>
-#f - funkcja f z rownania
+#f - funkcja f z rownania<br>
-#x0-wartosc poczatkowa
+#x0-wartosc poczatkowa<br>
-#t0-czas poczatkowy
+#t0-czas poczatkowy<br>
-#tk-czas koncowy
+#tk-czas koncowy<br>
-#h-krok czasowy
+#h-krok czasowy<br>
+<br>
-def EE(f,x0,t0,tk,h):
+def EE(f,x0,t0,tk,h):<br>
-     #generujemy wektor czasow
+     #generujemy wektor czasow<br>
-     t=np.arange(t0,tk,h)
+     t=np.arange(t0,tk,h)<br>
+<br>
-     #liczba krokow czasowych
+     #liczba krokow czasowych<br>
-     N=len(t)
+     N=len(t)<br>
+<br>
-     #wektor wynikowy
+     #wektor wynikowy<br>
-     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0))) # gdy mamy do czynienia w równaniem wektorowym
+     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0))) # gdy mamy do czynienia w równaniem wektorowym<br>
-     else: x=np.zeros(N) #dla przypadku skalarnego
+     else: x=np.zeros(N) #dla przypadku skalarnego<br>
+<br>
-     #wpisujemy wartosc poczatkowa
+     #wpisujemy wartosc poczatkowa<br>
-     x[0]=np.array(x0)
+     x[0]=np.array(x0)<br>
-     #index
+     #index<br>
-     i=1
+     i=1<br>
-     #petla glowna
+     #petla glowna<br>
-     while (i<N):
+     while (i<N):<br>
-         x[i]=np.array(x[i-1]+h*f(t[i-1],x[i-1]))
+         x[i]=np.array(x[i-1]+h*f(t[i-1],x[i-1]))<br>
-         i+=1
+         i+=1<br>
-     return t,x
+     return t,x<br>
 </syntaxhighlight>
 <!--</source>-->

@@ Linia 261: / Linia 261: @@
 [[Plik:oscy1.png]]
-Widać, że początkowo układ dochodzi do stanu regularnych oscylacji. Do analizy 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
+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">

@@ Linia 302: / Linia 302: @@
 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===
-Zacznijmu od implementacji podanej w treści zadania metody całkowania RK2
+Zacznijmy od implementacji podanej w treści zadania metody całkowania RK2
 <source lang="python">
 import numpy as np
@@ Linia 374: / Linia 374: @@
 [[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
+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])
@@ Linia 396: / Linia 396: @@
 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
+Okres trajektorii periodycznej możemy znaleźć na przykład w ten sposób
 <source lang="python">
 prog=140

@@ Linia 269: / Linia 269: @@
 </source>
-Interesować nas będzie amplituda drgań w funkcji częstość W. Wygenerujmy listę wartości W dla których będziemy liczyć amplitudę.
+Interesować nas będzie amplituda drgań w funkcji częstość W. Wygenerujmy listę wartości W, dla których będziemy liczyć amplitudę.
 <source lang="python">
@@ Linia 278: / Linia 278: @@
 amp=[amplituda(RK4(oscylator,[0.0,1.0],0.0,400.0,0.1)[1]) for oscylator.W in Omegas]
 </source>
-Wynik koncowy wyglada następująco
+Wynik końcowy wygląda następująco:
 <source lang="python">
 py.plot(Omegas,amp)

@@ Linia 239: / Linia 239: @@
 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 pierszego stopnia i zapisania go w postaci funkcji
+Zacznijmy od sprowadzenia równania drugiego stopnia do równania pierwszego stopnia i zapisania go w postaci funkcji
 <source lang="python">
@@ Linia 253: / Linia 253: @@
 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.
+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))
@@ Linia 261: / Linia 261: @@
 [[Plik:oscy1.png]]
-Widać że początkowo układ dochodzi do stanu regularnych oscylacji. Do analizy amplitudy interesuje nas jedynie koncowa 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
+Widać, że początkowo układ dochodzi do stanu regularnych oscylacji. Do analizy 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">

@@ Linia 192: / Linia 192: @@
 ==Runge-Kutta czwartego rzędu==
-Ostatnią metodą, którą omówimy jest najbardziej popularna metoda zwana w skrócie RK4. Metoda to 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
+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>,
@@ Linia 229: / Linia 229: @@
      return t,x
 </source>
 ==Przykłady==
 ===Zadanie - Wahadło matematyczne z tłumieniem i siłą wymuszającą===

@@ Linia 103: / Linia 103: @@
 ==Metoda Żabiego Skoku==
-W poprzedniej metodzie liczyliśmy wartość funkcji f w chwili <math>t_n</math>, która była pochodzną 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ą pochodzną 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>
+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>.
@@ Linia 109: / Linia 109: @@
 <math> k_2=f(t_n+h/2,x_n+h/2*k_1) </math>.
-Nastepnie używamy pochodznej <math> k_2 </math> zamiast pochodznej <math> k_1 </math> do obliczenia wartości funkcji w kolejnym kroku czasowym.
+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 nastepująco
+Przykładowa implementacja tej metody wygląda następująco
 <source lang="python">
 def leapfrog(f,x0,t0,tk,h):
@@ Linia 141: / Linia 141: @@
      return t,x
 </source>
-Rozwiązanie równania oscylatora tą metodą dla identycznych jak poprzednio czasów da nastepujące wyniki
+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)

@@ Linia 25: / Linia 25: @@
 ==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
+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łacjąc uzyskujemy
+przekształcając uzyskujemy
 <math> x(t+h)  \approx x(t)+ h \frac{d x(t)}{dt} </math>
@@ Linia 58: / Linia 58: @@
 def EE(f,x0,t0,tk,h):
-     #generujemy wektor czasow
+     #generujemy wektor czasowy
      t=np.arange(t0,tk,h)
@@ Linia 65: / Linia 65: @@
      #wektor wynikowy
-     if hasattr(x0, "__len__"): x=np.zeros((N,len(x0))) # gdy mamy doczynienia w równaniem wektorowym
+     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
@@ Linia 79: / Linia 79: @@
 </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
+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):
@@ Linia 86: / Linia 86: @@
      return np.array([xdot,-x])
 </source>
-Rozwiążmy to równanie z warunkeim początkowym [1.0,1.0] i od czasu od 0 do 100.
+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)
@@ Linia 92: / Linia 92: @@
 py.show()
 </source>
-rozwiązanie wygląda wówczas nastepująco
+rozwiązanie wygląda wówczas następująco
 [[Plik:img01.png]]
-Jeżeli zaś wydłuzymy czas symulacji to 1000 otrzymamy
+Jeżeli zaś wydłużymy czas symulacji do 1000 otrzymamy
 [[Plik:img2.png]]
-Amplituta 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łan 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
+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==