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

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"