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] 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].

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

"Programowanie dla Fizyków Medycznych"