# Runge Kutta módszerek (Explicit)
Készítette: Sykora Henrik - BME-GPK Műszaki Mechanikai Tanszék (sykora@mm.bme.hu)

Forrás: https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

### Kezdetiérték feladat (KDE):
$\dot{y}(t) = f(y,t)$  

### Runge Kutta módszerek:  
Közelítsük a megoldást diszkrét időpontokban ($t_n = nh$):  
$y_n \approx y(t_n)$  

$\begin{align}
y_{n+1} = y_{n} + h \sum_{i=1}^s b_i k_i,
\end{align}$  
ahol  
$\begin{align}
k_1 &= \left(t_n,y_n\right), \\
k_2 &= \left(t_n+c_2 h,y_n+\left(a_{21}k_1\right) h\right), \\
k_3 &= \left(t_n+c_3 h,y_n+\left(a_{31}k_1+a_{32}k_2\right)h\right), \\
&\;\vdots \\
k_s &= \left(t_n+c_s h,y_n+\left(a_{s1}k_1+\ldots+a_{s,s-1}k_{s-1}\right)h\right). 
\end{align}$  

Az $a_i$, $b_i$ és $c_i$ paramétereket az úgynevezett *Butcher tábla* (*Butcher tableau*) tartalmazza:

$\begin{equation}
\renewcommand\arraystretch{1.2}
\begin{array}
{c|cccc}
0\\
c_2 & a_{2,1}\\
c_3 & a_{3,1} & a_{3,2} \\
\vdots & \vdots& & \ddots \\
c_s & a_{s,1} & a_{s,2} & \dots & a_{s,1}\\
\hline
& b_1 & b_2 & \dots & b_{s-1} & b_s
\end{array}
\end{equation}$

#### Egyszerű esetek:
Euler módszer:  
$\begin{equation}
\renewcommand\arraystretch{1.2}
\begin{array}
{c|cccc}
0\\
\hline
& 1
\end{array}
\end{equation}$  
$\begin{align}
y_{n+1} &= y_{n} + h  \cdot 1 \cdot k_1,\\
k_1 &= f(t_n, y_n).
\end{align}$

Negyedrendű Runge-Kutta módszer:  
$\begin{equation}
\renewcommand\arraystretch{1.2}
\begin{array}
{c|cccc}
0\\
\frac{1}{2} & \frac{1}{2}\\
\frac{1}{2} &0 &\frac{1}{2} \\
1& 0& 0& 1\\
\hline
& \frac{1}{6} &\frac{1}{3} &\frac{1}{3} &\frac{1}{6} 
\end{array}
\end{equation}$  
$\begin{align}
y_{n+1} &= y_{n} + \frac{h}{6} \left(k_1+2k_2 + 2k_3 + k_4\right),\\
k_1 &= \left(t_n,y_n \right), \\
k_2 &= \left(t_n+\frac{1}{2} h, y_n+\frac{1}{2}k_1 h\right), \\
k_3 &= \left(t_n+\frac{1}{2} h, y_n+\frac{1}{2}k_2 h\right), \\
k_4 &= \left(t_n+ h, y_n+k_3 h\right). 
\end{align}$  

#### Megjegyzés: Adaptív Runge-Kutta módszerek
$\begin{equation}
\renewcommand\arraystretch{1.2}
\begin{array}
{c|cccc}
0\\
c_2 & a_{2,1}\\
c_3 & a_{3,1} & a_{3,2} \\
\vdots & \vdots& & \ddots \\
c_s & a_{s,1} & a_{s,2} & \dots & a_{s,1}\\
\hline
& b_1 & b_2 & \dots & b_{s-1} & b_s\\
& b_1^\ast & b_2^\ast & \dots & b_{s-1}^\ast & b_s^\ast
\end{array}
\end{equation}$  
Hibafüggvény ami alapján eldöntjük, hogy kell-e felezni/duplázni $h$ értékét:  
$e_{n+1} = y_n+1-y^\ast_n+1 = h \sum_{i=s}^s\left(b_i-b_i^\ast\right)k_i$

# Implementáció


```python
# Számításokhoz és az eredmények megjelenítéséhez szükséges csomagok behívása
import numpy as np
import matplotlib.pyplot as mpl
# Segédfüggvény
def getzero(ll):
  return np.zeros(len(ll))

# Runge Kutta lépés osztály
class RKstep:
  def __init__(self,a,b,c,h):
    self.a = a
    self.b = b
    self.c = c
    self.h = h
    self.nk = len(c)

  def dy(self,f,t,y): # f(t_n,y_n) -> k
      ks=[getzero(y) for i in range(self.nk)]
      for (i,ci) in enumerate(self.c):
        ks[i]= f(t+ci*self.h,y + self.h*sum(self.a[i][j]*ks[j] for j in range(i)))
      return self.h*sum(k*self.b[i] for (i,k) in enumerate(ks))

# RK4 módszer
def RK4(h):
  a=[[],[1/2],[0,1/2],[0,0,1/2]]
  b=[1/6,1/3,1/3,1/6]
  c=[0,1/2,1/2,1]
  return RKstep(a,b,c,h)

# Euler módszer
def Euler(h):
  a=[[]]
  b=[1.]
  c=[0.]
  return RKstep(a,b,c,h)


# Közönséges differenciálegyenlet kezdetiérték problémájának megadása/megoldásfüggvénye osztályként
class ODEProblem:
  def __init__(self,f,y0,tspan):
    self.f = f
    self.y0 = y0 if type(y0)==np.ndarray else np.array(y0)
    self.tspan = tspan
    
  def solve(self,method):
    t0,T = self.tspan
    ts = np.arange(t0,T+100*np.finfo(float).eps,method.h)
    ys=[self.y0 for (i,t) in enumerate(ts)]
    for i in range(len(ys)-1):
      ys[i+1] = ys[i] + method.dy(self.f,ts[i],ys[i])
    return (ts,ys)
```

# Alkalmazás

Pl.:  
$\dot{y}(t) = y(t), \quad y(0) = y_0, \quad t \in \left[0,T\right]$


```python
# Függvény és ODEProblem megadása
foo = lambda t,y: y#*np.sin(4*t)
y0=1.
T=3.
prob=ODEProblem(foo,[y0],(0.,T));

# Lépésköz + numerikus megoldás
h=0.1
# Runge Kutta megoldás
ts,ys=prob.solve(RK4(h))
#Euler módszer megoldás
tsE,ysE=prob.solve(Euler(h))

# Ábrázolás
mpl.plot(ts,np.exp(ts))
mpl.scatter(ts,[y[0] for y in ys],color='r',s=10)
mpl.scatter(tsE,[y[0] for y in ysE],color='g',s=10)
mpl.legend(['Pontos','RK4','Euler'])
mpl.show()
```

db949f461dad0d81c1f67e2ce15ec0542ef5a22e
ipynb
Jupyter Notebook
Python_RungeKutta.ipynb
HTSykora/Python_Bevezeto
c1d7dde5001d6218d0b4f6e53f304fb928fa92db
MIT
2021-03-01T17:31:09.000Z
2021-03-01T17:31:09.000Z
Python_RungeKutta.ipynb
HTSykora/Python_Bevezeto
c1d7dde5001d6218d0b4f6e53f304fb928fa92db
MIT
Python_RungeKutta.ipynb
HTSykora/Python_Bevezeto
c1d7dde5001d6218d0b4f6e53f304fb928fa92db
MIT
Qwen/Qwen-72B
1. YES
2. YES