.. _stap:

Jednodimenzioni direktni problem
==================================

Tanak štap od homogenog materijala je okružen izolacijom, tako da se promene temperature u štapu dešavaju samo kao posledica razmene toplote ka krajevima štapa i provođenja toplote duž štapa. Štap je jedinične dužine. Oba kraja su izložena mešavini vode i leda temperature 0. Početna temperatura na rastojanju :math:`x` od levog kraja štapa je :math:`sin{(\pi x)}`, kao što se vidi na :numref:`fdheat`.

.. _fdheat:

.. figure:: fdheat.png
    :width: 80%

    Eksperimentalna postavka problema provođenja toplote duž štapa. Na krajevima štapa nalazi se mešavina vode i leda. Štap je izolovan od uticaja spoljašnje sredine.

Radi poređenja, pokazaćemo sada kako se ovaj relativno jednostavan problem formuliše pomoću klasične metode konačnih razlika, a zatim ćemo rešiti direktni i inverzni problem opisane postavke NMPFZ metodom. 

Rešavanje metodom konačnih razlika (MKR)
-----------------------------------------

Parcijalna diferencijalna jednačina :math:numref:`eq:toplota1` modeluje temperaturu u bilo kojoj tački štapa u bilo kom vremenskom trenutku prema :cite:t:`recktenwald2004finite`. Ova jednačina se rešava metodom konačnih razlika, koja daje aproksimaciju rešenja za raspored temperature, primenjujući prostornu i vremensku diskretizaciju. Programska implementacija rešenja čuva temperaturu svake tačke diskretizacije u dvodimenzionoj matrici. Svaki red sadrži temperaturnu distribuciju štapa u nekom trenutku vremena. Štap je podeljen na :math:`n` delova dužine :math:`h`, pa stoga svaki red ima :math:`n+1` elemenata. Načelno, što je veće :math:`n`, manja je greška aproksimacije. Vreme od 0 do :math:`T` je podeljeno u :math:`m` diskretnih intervala dužine :math:`k`, pa stoga matrica ima :math:`m+1` redova, kao što je prikazano na :numref:`fdheat2`.

.. _fdheat1:

.. figure:: fdheat1.png
    :width: 80%

    Kako vreme teče, štap se hladi. Metoda konačnih razlika omogućava izračunavanje temperature u fiksnom broju tačaka u ravnomernim vremenskim intervalima. Smanjenje prostornog i vremenskog koraka uglavnom dovodi do preciznijeg rešenja.

.. _fdheat2:

.. figure:: fdheat2.png
    :width: 60%

    Diskretizacija jednačine provođenja toplote metodom konačnih razlika

Svaka tačka :math:`u_{i,j}` predstavlja element matrice koji sadrži temperaturu na poziciji :math:`i \cdot h`, u trenutku :math:`j \cdot k`. Na krajevima štapa je temperatura uvek nula. U početnom trenutku, temperatura u tački :math:`x` je, kao što je već rečeno, :math:`sin{(\pi x)}`. Algoritam ide korak po korak kroz vreme, koristi vrednosti iz trenutka :math:`j` da bi izračunao vrednosti u trenutku :math:`j+1`. Formula koja reprezentuje varijantu aproksimacije FTCS (*Forward Time Centered Space*) kao u :cite:t:`recktenwald2004finite` se ovde daje bez izvođenja i glasi:

.. math::
    :label: eq:diskretna

    u_{i,j+1} = R \cdot u_{i-1,j}+(1-2R) \cdot u_{i,j} + R \cdot u_{i+1,j},

gde je

.. math::

    R = \alpha \frac{k}{h^2}.

Celokupna analiza različitih eksplicitnih i implicitnih metoda data je u pomenutoj referenci, a ključni deo koda u programskom jeziku Pajton implementiran je na sledeći način:

.. code-block:: python
   :linenos:

    def heatFTCS(nt=10, nx=20, alpha=0.3, L=1, tmax=0.1):
        h = L / (nx - 1)
        k = tmax / (nt - 1)
        r = alpha * k / h**2

        x = np.linspace(0, L, nx)
        t = np.linspace(0, tmax, nt)
        U = np.zeros((nx, nt))

        # Početni uslov
        U[:, 0] = np.sin(np.pi * x / L)

        # Glavna petlja za MKR
        for m in range(1, nt):
            for i in range(1, nx-1):
                U[i, m] = r * U[i - 1, m - 1] + (1-2*r) * U[i, m-1] + r * U[i+1, m-1]

        # Egzaktno rešenje za poređenje
        ue = np.sin(np.pi * x / L) * \
            np.exp(-t[nt - 1] * alpha * (np.pi / L) * (np.pi / L))

Kao što detaljno objašnjava :cite:t:`recktenwald2004finite`, ako se MKR petlja formuliše eksplicitno kao što je to slučaj kod FTCS tehnike, mora se pažljivo izabrati vremenski i prostorni korak, kako bi numerička propagacija bila "brža" od fizičke. Rešenje koje se dobija pomoću MKR šeme se može videti na :numref:`ftcs`.

.. _ftcs:

.. figure:: ftcs.png
    :width: 60%

    Rešenje koje se dobija MKR metodom koristeći eksplicitnu FTCS tehniku u trenutku :math:`t=0.1s`

Ovaj problem ima i analitičko rešenje, pa je pogodan za testiranje različitih numeričkih metoda. To rešenje glasi:

.. math:: 
    :label: eq:analiticko1

    u(x,t) = \sin\left(\frac{\pi x}{L}\right) \cdot e^{-\frac{\alpha\pi^2}{L^2}t}.

ili u našem slučaju, kada je :math:`L=1`:

.. math:: 

    u(x,t) = \sin(\pi x) \cdot e^{-\alpha \pi^2 t}.

Eksplicitne tehnike poput FTCS ne garantuju konzistentnost rešenja koju garantuju implicitne tehnike kao što je BTCS (*Backward Time Centered Space*). MKR je ustaljeni pristup koji za većinu pravilno definisanih prostornih domena radi veoma dobro. Za ovako jednostavnu postavku kao što je jednodimenziono provođenje toplote i kada su svi parametri problema poznati (ovde je to :math:`\alpha`), MKR je često optimalna metoda. Međutim, kod većine problema iz prakse to nije slučaj. Hajde da razmotrimo kako da ovaj problem rešimo koristeći NMPFZ i direktno uporedimo sa MKR. 

Rešavanje pomoću NMPFZ
------------------------

Rešenje jednačine :numref:`fdheat` sa već postavljenim graničnim i početnim uslovima:

.. math:: 
    :label: eq:granicni1

    u(x=0,t)=u(x=1,t)=0, \, \forall t \\
    u(x,t=0) = sin{(\pi x)}

potražićemo pomoću NMPFZ pristupa. Iako je moguće da metode implementiramo direktno kao :cite:t:`raissi2019physics` koristeći okvir za duboko učenje kao što je `Tensorflow <https://www.tensorflow.org/>`_, ipak ćemo iskoristiti pomoć biblioteka koje dodatno apstrahuju NMPFZ entitete i omogućavaju korisniku da se fokusira na problem koji rešava. Ovaj primer rešićemo koristeći biblioteku `SCIANN <https://www.sciann.com/>`_ autora :cite:t:`haghighat2021sciann`. Postupak rešavanja objasnićemo direktno kroz programski kod:

.. code-block:: python
   :caption: NMPFZ - provođenje toplote
   :linenos:

    import numpy as np
    import sciann as sn
    from sciann.utils.math import diff, sign, sin, sqrt, exp
    from numpy import pi

    x = sn.Variable('x')
    t = sn.Variable('t')
    u = sn.Functional('u', [x,t], 3*[20], 'tanh')
    alpha = 0.3

    L1 = diff(u, t) - alpha * diff(u, x, order=2)

    TOLX = 0.011
    TOLT = 0.0011
    C1 = (1-sign(t - TOLT)) * (u - sin(pi*x))
    C2 = (1-sign(x - (0+TOLX))) * (u)
    C3 = (1+sign(x - (1-TOLX))) * (u)

    m = sn.SciModel([x, t], [L1, C1, C2, C3], 'mse', 'Adam')

    x_data, t_data = np.meshgrid(
        np.linspace(0, 1, 101), 
        np.linspace(0, 0.1, 101)
    )

    h = m.train([x_data, t_data], 4*['zero'], learning_rate=0.002, batch_size=256, epochs=500)
    
    # Test
    nx, nt = 20, 10
    x_test, t_test = np.meshgrid(
        np.linspace(0.01, 0.99, nx+1), 
        np.linspace(0.01, 0.1, nt+1)
    )
    u_pred = u.eval(m, [x_test, t_test])

Varijable :math:`x` i :math:`t` se na početku definišu na propisani način. Osnovni pojam koji se koristi u SCIANN biblioteci za apstrakciju NMPFZ je funkcional, koji je ovde označen sa :math:`u`, koji kao ulaz uzima :math:`x` i :math:`t`, ima 3 skrivena sloja sa po 20 neurona i kao aktivaciju svih tih neurona uzima funkciju hiperboličkog tangensa. Prvi sabirak kompozitne funkcije gubitka proizilazi iz same diferencijalne jednačine :math:numref:`eq:toplota1`. Kao što se vidi, za diferenciranje se koristi specijalni operator ``diff()`` iz biblioteke:

.. code-block:: python
    
    L1 = diff(u, t) - alpha * diff(u, x, order=2)

Najzanimljiviji i ne baš tako očigledan je način definisanja početnog uslova ``C1`` i graničnih uslova ``C2`` i ``C3``:

.. code-block:: python

    C1 = (1-sign(t - TOLT)) * (u - sin(pi*x))
    C2 = (1-sign(x - (0+TOLX))) * (u)
    C3 = (1+sign(x - (1-TOLX))) * (u)

Ovde je ``C1`` jednak nuli u svim tačkama uzorkovanja osim za :math:`t \le TOLT`. Tolerancije *TOLX* i *TOLT* su postavljene tako da "hvataju" prvu/poslednju vrstu ili kolonu kolokacionih tačaka, u zavisnosti šta je potrebno. Umesto funkcije znaka ``sign()``, mogu se koristiti i glatkije funkcije, kao što je hiperbolički tangens. NMPFZ model se formira pomoću ``SciModel`` konstruktora koji definiše i tip funkcije gubitka i algoritam optimizacije, tj. obučavanja:

.. code-block:: python

    m = sn.SciModel([x, t], [L1, C1, C2, C3], 'mse', 'Adam')

Obučavanje modela se pokreće metodom ``train()``, pri čemu se navode sledeći parametri:

#. **Skup kolokacionih tačaka za treniranje**.  Ovde je to pravilna ekvidistantna mreža tačaka po obe varijable. Horizontala je prostor, a vertikala vreme.
#. **Početne vrednosti komponenti funkcije gubitka**. Uobičajeno je da se na početku postave na nule.
#. **Stopa učenja**,
#. **Veličina batch-a**. Treba imati na umu da ako je broj tačaka domena u kojima se nameću granični uslovi značajno manji u odnosu na ukupan broj kolokacionih tačaka, parametar  optimizacije ``batch_size`` treba da bude podešen na veliki broj koji garantuje doslednu optimizaciju mini *batch*-a. U suprotnom, može da se desi da neki mini *batch*-evi ne dobiju nijednu kolokacionu tačku koja pripada graničnim uslovima i stoga ne generišu tačan gradijent za ažuriranje metodom gradijentnog spusta. 
#. **Broj epoha**.

.. _loss1:

.. figure:: loss1.png
    :width: 80%

    Istorija obuke jednodimenzionog modela provođenja toplote.

Tok obuke možemo da ispratimo kroz standardne *Tensorflow* objekte, kao što je ``h.history['loss']``, kao što se vidi na :numref:`loss1`. Pošto se završi obuka NMPFZ-a, možemo formirati testni skup tačaka slično kao što smo to učinili i sa kolokacionim tačkama i proveriti rezultate predikcije pozivom metode ``eval()`` na objektu istreniranog modela. Rezultat polja temperature duž štapa u trenutku :math:`t=0,1` i njegovo poređenje sa analitičkim rešenjem vidi se na :numref:`heat1`.

.. _heat1:

.. figure:: heat-pinn1.png
    :width: 80%

    Polje temperature duž štapa u trenutku :math:`t=0,1` dobijeno metodom NMPFZ.

Čisto praktično gledano, **NMPFZ rešenje jednostavnog direktnog problema kao što je ovaj i ne pruža nikakve posebne prednosti u odnosu na klasičnu MKR metodu**. Prvo, rešavanje duže traje i zahteva upošljavanje više računarskih resursa i zavisnosti u vidu dodatnih biblioteka za tenzorski račun. Dalje, specifikacija početnih i graničnih uslova kod NMPFZ ima svoje specifičnosti. Treće, neophodno je metodom probe i greške podesiti hiper-parametre modela, kao što su: broj skrivenih slojeva, broj neurona po sloju, aktivaciona funkcija, brzina učenja (*Learning Rate*) itd. Od izbora hiper-parametara konvergencija rešenja može značajno da zavisi. 

Sa druge strane, za razliku od MKR i MKE (*Metoda Konačnih Elemenata*), NMPFZ nam dozvoljava da problem definišemo čistim diferencijalnim jednačinama i proizvoljnim graničnim uslovima (Dirihleovi, Nojmanovi, periodični, skup tačaka). **Nema potrebe za specificiranjem algebarske veze između čvorova** (tj. kolokacionih tačaka u NMPFZ) i rešavanjem tako postavljenog sistema jednačina. Zahvaljujući ovoj činjenici, bilo koja nova fizika u vidu novog graničnog uslova ili promena u samoj diferencijalnoj jednačini može da se izvede veoma lako, omogućavajući brzu proveru hipoteza i izradu prototipova. 

Drugo, dok sve klasične metode proračun moraju da izvedu kroz vremenske korake (*time stepping*), **NMPFZ omogućava brzu inferenciju** na već obučenoj mreži za bilo koji vremenski trenutak :math:`t` postavljen na ulazu mreže. Za neke primene u realnom vremenu gde je brzina od ključnog značaja, ovo može da bude presudno. 

Treće, NMPFZ metodološki ne razlikuje **direktne probleme** (u kojima se rešava poznata diferencijalna jednačina) od **inverznih problema**, kod kojih su neki od parametara nepoznati, ali postoje dodatni uslovi iz kojih se nepoznati parametri mogu dobiti procesom treninga. U narednoj temi :ref:`stap_inverzni` demonstriraćemo jedan takav problem.