Matematikai algoritmusok és felfedezések 2.



Vendégelőadás:
Differenciálegyenletek numerikus módszerei és alkalmazásai



Fekete Imre
Alkalmazott Analízis és Számításmatematikai Tanszék
Matematikai Intézet

Hogyan modellezhetünk valós folyamatokat?


  • Principle based (klasszikus irányzat --> társtudományok elméletén alapuló)

  • Data-driven (manapság új és divatos irányzat --> adat, gépi tanulás / neurális háló alapú)

  • Hibrid (a fenti kettő kombinálása)

Milyen lehetőségek állnak rendelkezésünkre a klasszikus irányzaton belül?


  • Közönséges differenciálegyenletek (klasszikus irányzat --> társtudományok elméletén alapuló)
    • Differenciálegyenletek kurzus: megoldhatósági tételek, tételek speciális egyenletek és rendszerek kapcsán
    • Numerikus módszerek kurzus: időbeli diszkretizációs módszerek: Runge-Kutta, többlépéses módszerek és társaik

Milyen lehetőségek állnak rendelkezésünkre a klasszikus irányzaton belül?


  • Parciális differenciálegyenletek (manapság új és divatos irányzat --> adat, gépi tanulás / neurális háló alapú)
    • Funkcionálanalízis, Parciális differenciálegyenletek: mélyebb analízis eszköztár segítségével elméleti vizsgálat különböző típusú egyenletekhez (elliptikus, parabolikus, hiperbolikus)
    • Elliptikus, Időfüggő parciális differenciálegyenletek numerikus megoldási módszerei: véges differencia / véges térfogat / végeselem módszerek, speciális időintegrátorok

A teljes megvalósítási folyamat nagyon komplex



Kezdetiérték-feladatok


Legyen $G\subset \mathbb{R}\times \mathbb{R}^d$ összefüggű nyílt halmaz, $(t_0,u_0)\in \mathbb{G}$ egy adott pont, $f:G\to \mathbb{R}^d$ egy folytonos leképezés. A

\begin{cases} \frac{du(\cdot)}{dt}= f(\cdot,u)&\\ u(t_0)=u_0, & \end{cases}

feladatot kezdetiérték-feladatnak (Cauchy-feldatnak) nevezzük.

Numerikus módszerek - Bevezető


Explicit Euler (1768)

A közelítéshez a derivált fogalmát és a $\color{blue}{\text{jobboldali differenciahányadost}}$ használhatjuk:

$$ u'(t)=\lim_{h\to 0}\frac{u(t+h)-u(t)}{h}=f(t,u(t)) $$


$$ u'(t) \approx \color{blue}{\frac{u(t+h)-u(t)}{h}}=f(t,u(t)) $$


Cél: Minden rácspontban a pontos és numerikus közel legyen egymáshoz: $u(t_n) \approx y(t_n)$





Ekkor a javasolt numerikus módszer:

\begin{cases} \color{blue}{y_{n+1}= y_n + h f(t_n,y_n), \quad n=0,1,\ldots} &\\ \color{blue}{y_0=u_0} & \end{cases}

Mennyire jó közelítés az explicit Euler módszer?



A pontos megoldás adott rácspontjaiban Taylor-sorfejtés után: $$ \frac{u(t_n)+h\cdot u'(t_n) + \frac{h^2}{2!}\cdot u''(t_n)+\mathcal{O}(h^3)-u(t_n)}{h}=f(t_n,u(t_n)) $$

Egyszerűsítés és összevonás után a hibatag az alábbi: $$ \frac{\color{blue}{h^{1}}}{2}\cdot u''(\xi) $$

Implicit Euler


$\color{red}{\text{Baloldali differenciahányadost}}$ használva a numerikus módszer:

\begin{cases} \color{blue}{y_{n+1}= y_n + h} \color{red}{f(t_{n+1},y_{n+1})}\color{blue}{ \quad n=0,1,\ldots} &\\ \color{blue}{y_0=u_0} & \end{cases}



Taylor-sorfejtés és egyszerűsítés után a hibatag: $$ -\frac{\color{red}{h^1}}{2}\cdot u''(\xi) $$

Dahlquist-féle tesztfeladat és A-stabilitás


  Germund Dahlquist (1925-2005)
  Nemzetiség: Svéd

  Amiről híres: Teszt egyenlet (A-stabilitás), Többlépéses módszerek stabilitás elmélete, G-stabilitás, ...



A numerikus világ egyik, hanem a legfontosabb díját róla nevezték el: SIAM Germund Dahlquist Prize

Dahlquist-féle tesztfeladat


\begin{cases} u'(t)= \lambda u(t),& \lambda\in\mathbb{C}\\ u(0)=u_0, & \end{cases}

A megoldás: $u(t)=u_0e^{\lambda t}$. A tesztfeladatra alkalmazott tetszőleges egylépéses módszer esetén

$$y_{n+1} = R(h\lambda)y_{n}.$$
  • $\color{blue}{\text{Explicit Euler:}\ R(h\lambda) = 1 + h\lambda}$
  • $\color{red}{\text{Implicit Euler:}\ R(h\lambda) = 1 / (1 + h\lambda)}$

A-stabilitás


A tesztfeladat esetén a megoldás korlátosságához: $$Re(\lambda)\leq 0$$
Az $R(h\lambda)$ stabilitási függvényre vonatkozóan a stabilitási tartomány: $$S = \{ h\lambda\in\mathbb{C}: |R(h\lambda)|\leq 1 \}$$

Ezeket összevetve az A-stabilitáshoz kell: $$S\supset \mathbb{C}^{-}$$

In [1]:
## A NODEpy egy speciális nemlineáris közönséges differenciálegyenlet megoldó csomag
from nodepy import rk
explicit_euler = rk.loadRKM('FE')
explicit_euler.plot_stability_region()
Out[1]:
In [2]:
from nodepy import rk
implicit_euler = rk.loadRKM('BE')
implicit_euler.plot_stability_region()

The stability region is unbounded
Out[2]:

Konklúzió és egy egyszerű alkalmazás

$\color{blue}{\text{Explicit Euler}}$

  • $\color{blue}{\text{Előny: Számításmatematikai értelemben olcsó}}$
  • $\color{blue}{\text{Hátrány: Kicsi stabilitási tartomány (nem A-stabil módszer)}}$

$\color{red}{\text{Implicit Euler}}$

  • $\color{red}{\text{Előny: Nagy stabilitási tartomány (A-stabil módszer)}}$
  • $\color{red}{\text{Hátrány: Költséges (implicit és nemlineáris egyenlet megoldás)}}$

Logisztikus modellek: babszár és az olasz autók


Logisztikus növekedés (Verhulst-féle modell, 1838)

\begin{cases} N'(t)= rN(t)\bigg(1-\dfrac{N(t)}{K}\bigg)&\\ N(0)=N_0, & \end{cases}

ahol

  • $r\in\mathbb{R}^+$ a növekedési ráta
  • $K\in\mathbb{R}^+$ a környezet eltartóképessége
  • $N(t)>0$ a populáció mérete a $t$ időpontban

Babszár növekedés paraméterei

  • $K=83.4772$
  • $r=0.977361$
  • $N_0=1.53008$

In [22]:
run expliciteuler.py
In [23]:
def f(t, u):
    K = 83.4772 
    r = 0.977361
    return r*u-r/K*u**2
In [25]:
h, t_n, y = expliciteuler(f,1.53008,0,10,200)

import matplotlib.pyplot as plt
plt.plot(t_n, y, 'mo')
plt.title('Explicit Euler a babszarnövekedésre')
plt.show()

Autók száma Olaszországban

  • $K=18.516\cdot 10^6$
  • $r=0.239318$
  • $N_0=0.194123\cdot 10^6$

In [6]:
def f(t, u):
    K = 18.516*10e+6
    r = 0.239318
    return r*u-r/K*u**2
In [7]:
h, t_n, y = expliciteuler(f,0.194123*10e+6,0,80,200)

import matplotlib.pyplot as plt
plt.plot(t_n, y)
plt.title('Explicit Euler az olasz autó eladásokra 1950 - 1980 között')
plt.axis([0, 30, 0, 18*10e+6])
plt.show()

Implicit Euler implementálási kérdései

Newton módszer a $\color{blue}{g(x)=0}$ nemlineáris egyenlet megoldására:

$$x_{n+1} = x_n - \frac{g(x_n)}{g'(x_n)}, \quad n=0,1,\ldots,\ \text{és}\ x_0\ \text{adott}$$


Tipikus leállási feltételek:

  • maximális iterációszám
  • abszolút hiba $||x_{n+1} -x_n||<\text{TOL}$

Implicit Eulerre alkalmazva: $$\color{blue}{g(y_{n+1})=}\color{red}{y_{n+1}-y_{n}-hf(t_{n+1},y_{n+1})}\color{blue}{=0}$$

In [26]:
run impliciteuler.py
In [27]:
h, t_n, y = impliciteuler(1.53008,0,10,200)

import matplotlib.pyplot as plt
plt.plot(t_n,y)
plt.title('Implicit Euler a babszarnövekedésre')
plt.show()

Runge-Kutta (RK) módszerek (1901)


  Carl Runge (1856-1927)
  Nemzetiség: Német
  Doktori témavezető: Karl Weierstrass
  Amiről híres: Runge-Kutta módszerek


  Martin Kutta (1867–1944)
  Nemzetiség: Német (Lengyelországban született)
  Amiről híres: Runge-Kutta módszerek, Kutta feltétel aerodinamikában

A leghíresebb módszer ötlete - Az RK4 módszer


Az RK4 módszer kompakt implementálható alakja

$$ \begin{align} k_1 &= f(t_n,y_n) \\ k_2 &= f\left(t_n+\frac{1}{2}h,y_n+\frac{1}{2}hk_1\right) \\ k_3 &=f\left(t_n+\frac{1}{2}h,y_n+\frac{1}{2}hk_2\right) \\ k_4 &=f\left(t_n+h,y_n+hk_3\right) \\ y_{n+1} &=y_n+h\left(\frac{1}{6}k_1+\frac{1}{3}k_2+\frac{1}{3}k_3+\frac{1}{6}k_4\right) \end{align} $$

Az RK4 módszer tulajdonságai

  • Módszer rendje: 4
  • Stabilitási függvénye: $R(h\lambda) = 1 + (h\lambda) + \frac{(h\lambda)^2}{2!} + \frac{(h\lambda)^3}{3!} \frac{(h\lambda)^4}{4!}$
  • A módszer nem A-stabil

A zsákmány-ragadozó modell



Tekintsük a Lotka-Volterra modell (1926) modell konkrét alakját

$$P'_1(t)=\beta_1P_1(t)-\alpha_1P_1(t)P_2(t)$$$$P'_2(t)=-\beta_2P_2(t)+\alpha_2P_1(t)P_2(t),$$

ahol

  • $P'_1(t)$ a zsákmány (prey), míg $P'_2(t)$ a ragadozók (predator) időbeli változását jelentik
In [28]:
run RK4sys.py
In [29]:
def f(t, y):
    y1, y2 = y
    return [beta1*y1-alpha1*y1*y2, -beta2*y2+alpha2*y1*y2]

beta1 = 0.6
beta2 = 0.6
alpha1 = 0.3 * 10 ** -4
alpha2 = 1.5 * 10 ** -5
In [30]:
h, t, y = RK4sys(f, [90000, 2000], 0, 50, 10000)

first = y[:,0]
second = y[:,1]
import matplotlib.pyplot as plt
fig = plt.figure()
l1, l2 = plt.plot(t, first, t, second)
fig.legend((l1, l2), ('Zsákmány (Nyúl / Hal)', 'Ragadozó (Róka / Ember)'),loc ="upper right")

plt.show()
In [13]:
h, t, y = RK4sys(f, [90000, 2000], 0, 50, 10000)

fig = plt.figure()
plt.plot(y[:,0], y[:,1])
Out[13]:
[<matplotlib.lines.Line2D at 0x7fb3f8f18c70>]

Légköri konvekció és a "pillangó hatás"


Az egyszerűsített Lorenz modell (1963):

\begin{cases} x'(t)=\sigma(y(t)-x(t))&\\ y'(t)=rx(t)-y(t)-x(t)z(t)&\\ z'(t)=x(t)y(t)-bz(t)& \end{cases}

ahol

  • $\sigma=10$
  • $r=28$
  • $b=8/3$
In [31]:
def f(t, y):
    y1, y2, y3 = y
    return [sigma*(y2-y1), r*y1-y2-y1*y3, y1*y2-b*y3]

sigma = 10
b = 8/3
r = 28
In [32]:
h, t, y = RK4sys(f, [0, 1, 0], 0, 100, 10000)
first = y[:,0]
third = y[:,2]
fig = plt.figure()
plt.plot(first, third)
plt.xlabel('x')
plt.ylabel('z')
plt.title('Lorenz attraktor, x-z sík')
plt.show()

Az RK módszerek kompakt implementálható alakja



$$ \begin{align} k_1 &= f\left(t_n+c_1h,y_n+h[a_{11}k_1+\ldots+a_{1s}k_s]\right) \\ k_2 &= f\left(t_n+c_2h,y_n+h[a_{21}k_1+\ldots+a_{2s}k_s]\right) \\ &\vdots \\ k_s &= f\left(t_n+c_sh,y_n+h[a_{s1}k_1+\ldots+a_{ss}k_s]\right) \\ y_{n+1} &=y_n+h\left(b_1k_1+\ldots+b_sk_s\right) \end{align} $$

Butcher tabló


  John Butcher (1933 - )
  Nemzetiség: Új-zélandi
  Amiről híres: Butcher tabló, Butcher-fa, Runge-Kutta klub

Runge-Kutta módszerek Butcher tablója



$$ \begin{array} {c|cccc} c_1 & a_{11} & a_{12} & \ldots & a_{1s}\\ c_2 & a_{21} & a_{22} & \ldots & a_{2s}\\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_s & a_{s1} & a_{s2} & \ldots & a_{ss}\\ \hline & b_1 & b_2 & \ldots & b_s \end{array}\quad \quad \begin{array} {c|c} c & A\\ \hline & b^T \end{array} $$

Az EE, IE és RK4 Butcher tablói



$$ \begin{array} {c|c} 0 & 0\\ \hline & 1 \end{array}\quad\quad\quad \begin{array} {c|c} 1 & 1\\ \hline & 1 \end{array}\quad\quad\quad \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} $$

Adaptív (változó lépésközű) Runge-Kutta módszerek


Lokális hibabecslés és lépésköz frissítés


$$ \begin{array} {c|c} c & A\\ \hline & b^T\\ & \hat{b}^T \end{array}\quad \Longrightarrow \quad r_n \quad \Longrightarrow \quad h_{n+1}=\left(\frac{\text{TOL}}{r_n}\right)^{1/q}h_n $$

Brusselator probléma

Tekintsük az alábbi két reakciós autokatalikus reakció elméleti modelljét:

\begin{cases} u_1'(t)=A+u_1(t)^2u_2(t)-(B+1)u_1(t)&\\ u_2'(t)=Bu_1(t)-u_1(t)^2u_2(t) \end{cases}

ahol $A$ és $B$ pozitív konstansok.

In [16]:
from scipy.integrate import solve_ivp
In [17]:
def f(t,y):
    A = 1
    B = 3
    y1, y2 = y
    return [A + y1 ** 2 * y2 - (B + 1) * y1, B * y1 - y1 ** 2 * y2]
In [18]:
sol = solve_ivp(f, [0, 20], [1.5, 3],rtol = 10 ** -10, atol = 10 ** -10)
import matplotlib.pyplot as plt
fig = plt.figure()
a1, a2 = plt.plot(sol.t, sol.y[0], sol.t, sol.y[1])
fig.legend((a1,a2), ('Első komponens','Második komponens'),'upper right')
plt.show()
In [19]:
import numpy as np

h = np.zeros(len(sol.t)-1)
print(len(sol.t))
for i in range (0,len(sol.t)-1):
    h[i] = sol.t[i+1]-sol.t[i]

print(np.amin(h[0:-1]))
print(np.amax(h))
plt.semilogy(sol.t[1:-1], h[0:-1], 'b*')

647
0.0037172656319675534
0.0998352440243071
Out[19]:
[<matplotlib.lines.Line2D at 0x7fb3f9a5dd30>]

Chemakzo feladat


Irányításelmélet


  Gustaf Söderlind (1952 – )
  Nemzetiség: Svéd
  Amiről híres: Adaptív RK módszerek, merevség, logaritmikus norma

Chemakzo feladat


$h_{n+1}=\left(\frac{\text{TOL}}{r_n}\right)^{\beta_1/q}\left(\frac{\text{TOL}}{r_{n-1}}\right)^{\beta_2/q}\left(\frac{h_n}{h_{n-1}}\right)^{-\alpha}h_n$


Szabályozó $\beta_1$ $\beta_2$ $\alpha$
I 1 0 0
PI3333 2/3 -1/3 0
H211PI 1/6 1/6 0
H211b 1/b 1/b 1/b

Közönséges differenciálegyenletek és gépi tanulás


Trajektória tanulásra vonatkozó szimulációs eredmények

Még több eszköztárral további izgalmas feladatok várhat ránk!

Köszönöm szépen a figyelmet!



ai.elte.hu