TI/Programowanie dla Fizyków Medycznych/RRZ
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]
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]),')'
