\documentclass[a4paper,fleqn]{article}
\usepackage[latin2]{inputenc}
\usepackage[magyar]{babel}
\usepackage{amssymb,graphicx,hhline}
\def\idr{\rm{id}_\mathbb{R}}
\def\d{\mbox{\rm d}}
\def\ii{{\rm i}}
\def\e{\mbox{\rm e}}
\def\be{\begin{equation}}
\def\ee{\end{equation}}

% Alap tav a lap bal illetve folso szeletol 2.54cm (1inch). Ezt valtoztatja a hoffset es a voffset.
% Az A4 lap szelessege 21cm, magassaga 29.7cm.
\setlength{\hoffset}{-40 pt}
\setlength{\voffset}{-30 pt}
\setlength{\topmargin}{0 pt}
\setlength{\textwidth}{400 pt}
\setlength{\textheight}{640 pt}

\begin{document}
\title{Differenciálegyenletek numerikus megoldása\\Jegyzet-féleség...}
\author{Csanád Máté}
\maketitle

\tableofcontents
\pagebreak
\section{Közönséges differenciálegyenletek}
\subsection{Diffegyenletek diszkretizációja, diszkretizációs lépés}
Egy $f$ függvény deriváltját az
\be f'(x)=\lim_{\Delta x \rightarrow 0}\frac{f(x+\Delta x)-f(x)}{\Delta x}\ee
formulával definiáljuk. Ezt a deriváltat közelítőleg számolhatjuk numerikusan, a probléma által megkövetelt pontosságnak
megfelelően kicsi $\Delta x$-et választva, a határérték elhagyásával:
\be f'(x)\approx\frac{f(x+\Delta x)-f(x)}{\Delta x}.\label{e:disq1}\ee
Ezt nevezzük a differenciálás diszkretizálásának. A fenti formula elsőrendben pontos, azaz a hiba arányos a lépésköz,
azaz $\Delta x$ négyzetével. Vannak más, pontosabb formulák, amelyekben az $f$ függvény $x$-en és $x+\Delta x$-en kívül
más pontokban is vett értékeinek segítségével számítjuk ki az $f'$ deriváltat az $x$ pontban, azaz $f'(x)$-et.
Általánosságban így néz ki egy diszkretizáció:
\be f'(x)\approx D(f(x_1),f(x_2),...,f(x_n);\Delta x),\ee
ahol $x_i$-k $x$ egy környezetében lévő értékek, ezeknek egy kombinációját jelöli $D$.

Amennyiben egy elsőrendű differenciál-egyenletet akarunk megoldani, többnyire az alábbi egyenletünk van:
\be h(f'(x),f(x),x)=0\label{e:deq1},\ee
például
\be f'(x)-\lambda f(x)=0\ee
vagy
\be f'(x)-\frac{2x}{x^2-1}=0.\ee
Az (\ref{e:deq1}) egyenletet a diszkretizáció segítségével a
\be h(D(f(x_1),f(x_2),...,f(x_n);\Delta x),f(x),x)=0\label{e:deq2}\ee
egyenletté alakíthatjuk. Ekkor vesszük az $\{x,x_1,...,x_n\}$ halmazból a legnagyobb értéket, legyen ez $x_k$, és erre
megoldjuk az egyenletet. Ekkor a következőre jutunk:
\be f(x_k)=H(x_1,..,x_{k-1},x_{k+1},..,x_n,x).\label{e:deq3}\ee

A fenti (\ref{e:deq3}) egyenlet segítségével \emph{léphetünk} egyet előre a megoldásban, ugyanis kisebb $x_i$ értékek
ismeretében kiszámítottuk egy következő $x_k$ helyen a függvény értékét. Azaz, ha ismerjük a függvény értékét megfelelő
mennyiségű kezdeti pontban (ezek a kezdeti feltételek), tetszőleges további helyeken kiszámíthatjuk a rá vonatkozó
diffegyenlet segíségével!

Lássunk most a fenti sémának megfelelő módszereket!


\subsection{Euler-módszer}
Az Euler módszer a fenti (\ref{e:disq1}) diszkretizációs formulán alapul. Ha van egy
\be f'(x)=H(f(x),x)\label{e:deq_euler}\ee
egyenletünk, akkor $f'(x)$-et meghatározzuk (\ref{e:disq1})-ből, és behelyettesítjük (\ref{e:deq_euler})-be.
Ekkor a következőre jutunk:
\be \frac{f(x+\Delta x)-f(x)}{\Delta x}=H(f(x),x).\ee
Innen pedig
\be f(x+\Delta x)=f(x)+\Delta x H(f(x),x)!\ee
Így $f(x_0)=f_0$ ismeretében kiszámíthatjuk a $f(x+n\Delta x)$ mennyiségek közelítő értékét.
Ezt az
\be f'(x)=-f(x) \ee
diffegyenlet esetére a következőképpen programozhatjuk be:

\bigskip
\noindent
\verbß  #include <iostream>ß\\
\verbß  #include <cmath>ß\\
\verbß  #include <fstream>ß\\
\verbß  ß\\
\verbß  using namespace std;ß\\
\verbß  ß\\
\verbß  int main()ß\\
\verbß  {ß\\
\verbß    double f0, f1;ß //a függvényérték az első és második pontban\\
\verbß    double dx, x;ß //a futó változó és a lépésköz\\
\verbß    int Nx;ß //a lépesek száma\\
\verbß  ß\\
\verbß    Nx = 1000;ß\\
\verbß    dx = 0.01;ß //így $1000\times0.01=10$-ig jutunk majd el\\
\verbß  ß\\
\verbß    f0 = 1;ß //a kezdeti érték megadása\\
\verbß    x = 0;ß //$x$-ben nulláról indulunk\\
\verbß    ofstream outfile;ß\\
\verbß    outfile.open("result.out",ios::out);ß //a kiemeneti fájl\\
\verbß                                        ß //megnyitása írásra\\
\verbß    for(int i=1;i<Nx;i++)ß\\
\verbß      {ß\\
\verbß        outfile << x << " " << f0 << endl;ß //írás a kimeneti fájlba\\
\verbß        x += dx;ß //a futó változó inkrementálása\\
\verbß        f1 = f0+dx*(-f0);ß //a diszkretizációs lépés\label{p:euler}\\
\verbß        f0 = f1;ß //az új függvényérték visszarása az előző helyére\\
\verbß      }ß\\
\verbß  ß\\
\verbß    return 0;ß\\
\verbß  }ß\\

A program lefordítása és lefuttatása után elindíthatjuk a \textsc{gnuplot} nevű programot, és ábrázolhatjuk az eredményt.
Ezt a következőképpen tehetjük meg:

\bigskip
\noindent
\verbß  csanad@itl84:~/temp> g++ program.cc -o program.exeß //Fordítás\\
\verbß  csanad@itl84:~/temp> ./program.exeß //Futtatás\\
\verbß  csanad@itl84:~/temp> gnuplotß //A \textsc{gnuplot} elindítása\\
\verbß  gnuplot> plot 'result.out' u 1:2 w lß //Ábrázolás vonalakkal\\

Az eredmény az \ref{f:res_euler}. ábrán látható.

\begin{figure}
\begin{center}
  \includegraphics[width=200pt,angle=270]{result.eps}\\
  \caption{Az Euler-módszer használatának az eredménye}\label{f:res_euler}
\end{center}
\end{figure}

\subsection{Runge-Kutta módszer}
A fenti Euler-módszer hibája a lépésköz mégyzetével arányos, azaz
\be \Delta f \sim \Delta x^2. \ee
A Runge-Kutta módszer alkalmazásával ennél kisebb hiba érhető el, avagy megfordítva, ugyanakkora hiba nagyobb lépésköz
alkalmazásával kapható, azaz ugyanolyan távolságba kevesebb lépéssel juthatunk el, értékes processzor-időt takarítva
meg ezzel. A módszer azon alapul, hogy a függvény deriváltját másik pontban tekintjük, mint ahonnan indulunk.

Továbbra is az (\ref{e:deq1}) diffegyenletet akarjuk megoldani. Ehhez bevezetjük a
\be \Delta f_1 = \Delta x h(x_n,f_n)\textnormal{ és } \ee
\be \Delta f_2 = \Delta x h(x_n+\frac{\Delta x}{2},f_n+\frac{\Delta f_1}{2}) \ee
mennyiségeket. Ezekkel pedig a
\be f_{n+1} = f_n + \Delta f_2 \ee
lépést tesszük meg. Ha utánaszámolunk a keresett $f$ függvény Taylor-sorának segítségével, látjuk, hogy ennek
a módszernek a hibája a lépésköz köbével ará\-nyos, azaz harmadrendben pontos. Egyes tanulmányokban,
programokban ezt a módszert továbbfejlesztett sokszögek módszerének hívják, és bevezetnek egy másikfajta
(a Heun féle) másodrendű Runge-Kutta módszert, ahol
\be \Delta f_1 = \Delta x h(x_n,f_n)\textnormal{ és } \ee
\be \Delta f_2 = \Delta x h(x_n+\Delta x,f_n+\Delta f_1) \ee
mennyiségeket. Ezekkel pedig a
\be f_{n+1} = f_n + \frac{\Delta f_1+\Delta f_2}{2} \ee
lépést teszik meg. Ez ugyanahoz a pontossághoz vezet.

\label{p:rk4}
A leggyakrabban a negyedrendű Runge-Kutta módszer használják, ez pedig
a következőképpen működik. Itt a következő négy segéd-mennyiséget vezetjük be:
\be \Delta f_1 = \Delta x h(x_n,f_n), \ee
\be \Delta f_2 = \Delta x h(x_n+\frac{\Delta x}{2},f_n+\frac{\Delta f_1}{2}), \ee
\be \Delta f_3 = \Delta x h(x_n+\frac{\Delta x}{2},f_n+\frac{\Delta f_2}{2})\textnormal{ és} \ee
\be \Delta f_3 = \Delta x h(x_n+\Delta x,f_n+\Delta f_3). \ee
Ezekkel pedig:
\be f_{n+1} = f_n + \frac{\Delta f_1+2\Delta f_2+2\Delta f_3+\Delta f_4}{6}.\ee

\subsubsection{A Runge-Kutta módszer a gyakorlatban}
Ha másodrendű egyenletünk van, azaz például az
\be \ddot I + a \dot I + b I = c \ee
egyenlet, akkor ebből elsőrendű egyenletrendszert csinálnunk, azaz
\be \dot I = V \ee
\be \dot V = c - a V - b I. \ee
Itt a Runge-Kutta lépéseket párhuzamosan kell elvégeznünk. A diszkrecitációs lépés programozva a követekzőképpen nézne
ki:

\bigskip
\noindent
\verbß  I0=0;ß //kezdeti\\
\verbß  V0=1;ß //értékek\\
\verbß  ß\\
\verbß  kI1=dt*V0;ß //első segéd-mennyiség \verbß  Iß-re\\
\verbß  kV1=dt*(c-a*V0-b*I0);ß //első segéd-mennyiség \verbß  Vß-re\\
\verbß  ß\\
\verbß  kI2=dt*(V0+kV1/2);ß //második segéd-mennyiség \verbß  Iß-re\\
\verbß  kV2=dt*(c-a*(V0+kV1/2-b*(I0+kI1)/2);ß//második segéd-mennyiség \verbß  Vß-re\\
\verbß  ß\\
\verbß  kI3=dt*(V0+kV2/2);ß\\
\verbß  kV3=dt*(c-a*(V0+kV2/2-b*(I0+kV2/2));ß\\
\verbß  ß\\
\verbß  kI4=dt*(V0+kV3);ß\\
\verbß  kV4=dt*(c-a*(V0+kV3/2-b*(I0+kI3));ß\\
\verbß  ß\\
\verbß  ßV1=V0+kV1/6+kV2/3+kV3/3+kV4/6;\\
\verbß  ßI1=I0+kI1/6+kI2/3+kI3/3+kI4/6;\\
\verbß  ß\\
\verbß  ßV0=V1; //az új függvényérték visszamásolása az előzőbe\\
\verbß  ßI0=I1;\\

Azaz a fenti sorokat lehetne beírni a \ref{p:euler}-ben szereplő diszkretizációs lépés helyére.

Ennek a módszernek a hibája arányos a lépésköz ötödik hatványával.

\subsection{Adaptív lépéshossz-változtatás}

Ha egy olyan diffegyenletünk van, amelynek a megoldása egyes intervallumokon gyorsan oszcillál, míg máshol kellemesen
lapos, akkor nehéz megválasztani a lépésközt, mivel ha nagyot választunk, az a gyorsan oszcilláló tartományon rossz
eredményt ad, ha viszont kis lépésközt választunk, lassan haladunk.

Erre a problémára ad egy lehetséges megoldást az adaptív lépéshossz-változ\-tatás, melynek működésének lényegét az
alábbiakban láthatjuk. Tegyük fel, hogy egy adott $f_n$-ből indulunk, $\Delta x$ lépésközzel! Azt szeretnénk, hogy
a hiba ne lépjen túl egy előre megadott $\widetilde{\Delta f}$ értéket. Ekkor a következőket tesszük:
\begin{enumerate}
\item Megteszünk egy lépést $\Delta x$ lépésközzel, így $y_n$-ből $y_{n+1}$-be jutunk.
\item Megteszünk egy lépést $\frac{\Delta x}{2}$ lépésközzel, így $y^{(1/2)}_{n+1}$-be jutunk.
\item $y^{(1/2)}_{n+1}$-ből még egy lépést megteszünk $\frac{\Delta x}{2}$ lépésközzel, így $y^{(1/2)}_{n+2}$-be jutunk.
\item Összehasonlítjuk $y_{n+1}$-t és $y^{(1/2)}_{n+2}$-t. Ez a két érték közelítőleg egyenlő kellene, hogy legyen.
\item A két fenti érték alapján kapunk egy becslést a hibára:
\be \Delta f = \left|y_{n+1}-y^{(1/2)}_{n+2}\right|. \ee
\item Mivel tudjuk, hogy $\Delta f \sim \Delta x^5$ negyedrendű Runge-Kutta módszer esetén, ebből a kívánt hibához
szükséges lépésköz:
\be \widetilde{\Delta x}=\Delta x \left|\frac{\widetilde{\Delta f}}{\Delta f}\right|^{\frac{1}{5}} \ee
\item Innen két lehetőség van:
\begin{enumerate}
\item $\widetilde{\Delta x} \geq \Delta x$. Ekkor a következő lépés már lehet a nagyobb $\widetilde{\Delta x}$ méretű
\item $\widetilde{\Delta x} < \Delta x$. Ekkor az éppen elvégzett lépést az új, $\widetilde{\Delta x}$ lépésközzel újra
el kell végezni, és a következőnél is ezt használni.
\end{enumerate}
\item Vissza az első lépésre, amíg el nem értük a szükséges $x$ távolságot.
\end{enumerate}
\label{p:adapt}
Ezt az eljárást úgy programozhatjuk be, hogy definiálunk egy ``lépő" függ\-vényt, és azzal lépegetünk. Ez a
következőképpen néz ki mondjuk az
\be f'(x)=-f(x)*x \ee
diffegyenlet esetére:\\
\verbß  double step(double x0, double f0, double dx)ß\\
\verbß  {ß\\
\verbß    k1=dx*(-f0*x0);ß\\
\verbß    k2=dx*(-(f0+k1/2)*(x0+dx/2));ß\\
\verbß    k3=dx*(-(f0+k2/2)*(x0+dx/2));ß\\
\verbß    k4=dx*(-(f0+k3)*(x0+dx));ß\\
\verbß    ß\\
\verbß    return f0+k1/6+k2/3+k3/3+k4/6;ß\\
\verbß  }ß\\
A \verb"main()" függvénybe pedig a következők kerülnek:
\verbß  ß\\
\verbß  int main()ß\\
\verbß  {ß\\
\verbß    double f0, f1, f11, f12, df, dx;       ß //változók deklarálása\\
\verbß    double df0=0.01;                       ß //az elvárt pontosság\\
\verbß    double dx =0.01;                       ß //a kezdeti lépésköz\\
\verbß    for(int i=0;i++;i<imax)ß\\
\verbß      {ß\\
\verbß         f1 = step(x0, f0, dx);            ß //egész lépés\\
\verbß         f11= step(x0, f0, dx/2);          ß //első fél lépés\\
\verbß         f12= step(x0+dx0/2, f0+f11, dx/2);ß //második fél lépés\\
\verbß         df = fabs(f1-f12);                ß //becslés a hibára\\
\verbß         dx0= dx*pow(fabs(df0/df),0.2);    ß //a \verb"df0" hibához\\
\verbß                                             ß szükséges lépésköz\\
\verbß         if(dx0<dx)ß\\
\verbß           f1 = step(x0, f0, dx);          ß //az előző lépés megismétlése,\\
\verbß                                             ß már a kisebb lépésközzel\\
\verbß         dx = dx0;ß\\
\verbß         f0 = f1;ß\\
\verbß         x0+= dx;ß\\
\verbß         outfile << x0 << " " << f0 << endl;ß//írás a kimeneti fájlba\\
\verbß      }ß\\
\verbß  }ß\\
Még hasznosabb, ha még a diffegyenlet jobb oldalán álló függvényt is külön definiáljuk, azaz bevezetjük a\\
\verbß  double h(double f, double x)ß\\
\verbß  {ß\\
\verbß    return -f*x;ß\\
\verbß  }ß\\
függvényt, és ekkor a \verb"step(double,double,double)" függvénybe az alábbiakat írhatjuk:\\
\verbß  k1=dx*h(f0,x0);ß\\
\verbß  k2=dx*h(f0+k1/2,x0+dx/2);ß\\
\verbß  k3=dx*h(f0+k2/2,x0+dx/2);ß\\
\verbß  k4=dx*h(f0+k3,x0+dx);ß\\
Sőt, lehet két lépést is definiálni, egy \verb"rk_step(double,double,double)"-t és egy \\
\verb"euler_step(double,double,double)"-t, ahol az előbbi a fentieknek megfelelően, az utóbbit pedig egyszerűen az
Euler-módszer alapján egy\\
\verbß    return f0+dx*h(f0,x0);ß\\
sorral definiálhatjuk.

\section{Parciális differenciálegyenletek}
\subsection{Diszkretizáció}
Egy egyszerű parciális diffegyenlet a következő például:
\be
\frac{\partial \rho}{\partial t} + v \frac{\partial \rho}{\partial x} = 0
\ee
Ebben kétféle derivált is szerepel, és mindkettőt diszkretizálnunk kell. Ezt megtehetjük az Euler-módszerrel.
Itt viszont két indexe lesz a diszkretizált változónak, egy az idő-rácson, egy pedig a tér-rácson vett helyzetét adja
meg. Ennek az egyenletnek egyetlen paramétere van, méghozzá $v\Delta t / \Delta x$, ezt elnevezzük majd $C$-nek, ez a
Courant-szám. A diszkretizát egyenlet a követekzőképpen fog kinézni:
\be
\rho_j^{n+1} - \rho_j^n + v \frac{\Delta t}{\Delta x} \left(\rho_{j+1}^n - \rho_j^n\right) = 0
\ee
Igen egyszerű megoldást kaptunk így, ugyanis kifejezzük $\rho_j^{n+1}$-et a fenti egyenletből:
\be
\rho_j^{n+1} = \rho_j^n - C \left(\rho_{j+1}^n - \rho_j^n\right) \label{e:pdisq}
\ee
Így $\forall i$-re ismert $\rho^n_i$ esetén kiszámíthatjuk $\rho^{n+1}_i$-t $\forall i$-re. Természetesen szükségünk van
valami peremfeltételre, hogy a térbeli határokon hogyan viselkedjen a diszkretizált változó, ugyanis amikor a határon
vagyunk, azaz $j=j_{max}$, a $j+1$-edig indexnek nincs értelme. Két tipikus választás van: vagy a szélső ponton fixálom
nullára, egyre vagy bármely értékre a megoldást, vagy peridódikus határfeltételt veszek, azaz a $j_{max+1}$-edik pontot
azonosítom az elsővel. Természetesen mindkét esetben a kezdeti feltételnek, azaz $\rho^0_i$-nek is teljesítenie kell
a peremfeltételt.

\subsection{Stabilitás}
\label{ss:stab}
Sajnos a parciális differenciálegyenleteknél azonban nem ilyen egyszerű az élet, a fenti diszkretizáció ritkán műdödik.
Ugyanis a ``valódi" megoldáshoz mindig hozzáadhatunk egy olyan megoldást, ami a rácspontokon nulla, de közbe oszcillál.
Ez a rész-megoldás minden határon túl felerősödhet, és elnyomhatja az eredeti megoldást. Ezért egy speciális
diszkretizáció csak a $\Delta x$, $\Delta t$ és az egyéb paraméterek speciális megválasztása esetén működik. Azt úgy
vizsgáljuk, hogy bevezetjük a
\be
\rho_j^n = G^n \e^{\ii k \Delta x j} \label{e:stab}
\ee
megoldást. Ennek hullámszáma $k$, amelynek értékétől függő a megoldás perió\-dicitása. $G$ a hullám amplitúdóját jelöli,
és ha ez (vagyis ennek abszolútértéke) nagyobb egynél, akkor a megoldás az időben előrehaladva minden határon túl nő.

Vizsgáljuk meg, az (\ref{e:pdisq}). egyenletben szereplő diszkretizációnak milyen fel\-tételekkel megoldása
(\ref{e:stab})! Ehhez egymásba helyettesítjük őket, egyszerűsítünk $G^n\exp^{\ii k \Delta x j}$-vel, és ami marad:
\be
G = 1-C(\e^{\ii k \Delta x}-1)
\ee
Az ennek az egyenletnek megfelelő $G$, $C$ és $\Delta x$ értékek esetén a fenti, (\ref{e:stab})-ban szereplő függvény
megoldás lesz. Most lássuk a $|G|<1$ feltételt!
\be
|G| = (1+C-C \cos (k \Delta x))^2 + C^2 \sin ^2 (k \Delta x)
\ee
Ezen láthatjuk, hogy $G$ abszolútértéke mindig nagyobb mint egy, azaz ez a diszkretizáció mindig instabil.
Most lássunk speciális egyenltekre speciális diszkretizációkat!

\subsection{Kontinuitási egyenlet}
A fent említett parciális differenciálegyenlet a kontinuitási egyenlet. Erre nézünk most speciális diszkretizációkat.
\subsubsection{FTCS}
Az FTCS (Forward in Time, Centered in Space) módszer a követezőképpen néz ki:
\be
\frac{\rho_j^{n+1} - \rho_j^n}{\Delta t} + v \frac{\rho_{j+1}^n - \rho_{j-1}^n}{2 \Delta x} = 0
\ee
Ezt az \ref{ss:stab}-ben említett módszerrel megvizsgálva a következőre jutunk:
\be
G = 1 - C \ii \sin (k \Delta x)
\ee
Tehát ez feltétel nélkül instabil.
\subsubsection{Upwind}
Az upwind módszerben a térbeli deriváltat nem centráljuk, hanem áramlásirányból vesszük, azaz:
\be
\frac{\rho_j^{n+1} - \rho_j^n}{\Delta t} + v \frac{\rho_{j}^n - \rho_{j-1}^n}{\Delta x} = 0
\ee
vagy
\be
\frac{\rho_j^{n+1} - \rho_j^n}{\Delta t} + v \frac{\rho_{j+1}^n - \rho_{j}^n}{\Delta x} = 0
\ee
ha $v<0$.
Ekkor a stabilitásra a következőt kapjuk:
\be
G = 1-|C| + |C| \e^{\ii k \Delta x}.
\ee
Ez a komplex szám az egységkörön belül van, ha $|C| < 1$, azaz ekkor stabil a diszkretizáció.
\subsubsection{Leapfrog}
\label{sss:lf}
Itt
\be
\frac{\rho_j^{n+1} - \rho_j^{n-1}}{2 \Delta t} + v \frac{\rho_{j+1}^n - \rho_{j-1}^n}{2 \Delta x} = 0.
\ee
Ez a módszer szintén $|C| < 1$ esetén stabil, viszont térben és időben is másodrendűen pontos, szemben az előzőekkel.
\subsubsection{Lax-Wendorff}
\label{sss:lw}
Itt a kontinuitási egyenlet kétszer alkalmazzuk:
\be
\frac{\partial^2 \rho}{\partial t^2} + v^2 \frac{\partial^2 \rho}{\partial x^2} = 0.\label{e:lw}
\ee
Mivel a $\rho$ Taylor-sorából az
\be
\frac{\rho_j^{n+1} - \rho_j^{n}}{\Delta t} = \frac{\partial \rho}{\partial t}
+ \frac{\Delta t}{2} v^2 \frac{\partial^2 \rho}{\partial x^2} + ...
\ee
következik, a (\ref{e:lw}) egyenletet kihasználva a következőt írjuk fel:
\be
\frac{\partial \rho}{\partial t} = \frac{\rho_j^{n+1} - \rho_j^{n}}{\Delta t}
- \frac{\Delta t}{2} v^2 \frac{\partial^2 \rho}{\partial x^2}.
\ee
Így a következő, pontosabb diszkretizációt kapjuk:
\be
\rho_j^{n+1} = \rho_j^{n} -  \frac{C}{2} (\rho_{j+1}^n - \rho_{j-1}^n)
+ \frac{C^2}{2}(\rho_{j+1}^n  - 2 \rho_j^{n} + \rho_{j-1}^n).
\ee
Ez a módszer térben harmadrendig pontos, és szintén $|C| < 1$ esetén stabil.
\subsubsection{McCormack}
Ennél a módszernél két segéd-változót vezetünk be:
\begin{eqnarray}
\rho_j^\#&=&\rho_j^n-C(\rho_j^n-\rho_{j-1}^n)\label{e:contmc1}\\
\rho_j^{\#\#}&=&\rho_j^n-C(\rho_{j+1}^n-\rho_j^n),\label{e:contmc2}
\end{eqnarray}
azaz egyszer jobbról, egyszer balról vesszük a térbeli deriváltat. Itt ismét bevezettük $C$-t, a Courant-számot.
Ezekkel aztán kifejezzük az $n+1$-edik időpontban vett $\rho$-t:
\be
\rho_j^{n+1} = \frac{\rho_j^\#+\rho_j^{\#\#}}{2}\label{e:contmc}
\ee
Ha az a (\ref{e:contmc1})-(\ref{e:contmc1}) egyenleteket behelyettesítjük (\ref{e:contmc})-ba, akkor a következőt kapjuk:
\be
\rho_j^{n+1} = \rho_j^{n} - \frac{C}{2} (\rho_{j+1}^n-\rho_{j-1}^n) +
\frac{C^2}{2} (\rho_{j+1}^n - 2\rho_j^{n} + \rho_{j-1}^n)
\ee
Ez a módszer is $|C|<1$ esetén stabil, viszont egyszerűen számolható. A kontinuitási egyenlet esetén ez a módszer
megegyezik az \ref{sss:lw}-ben tárgyalt Lax-Wendorff módszerrel.
\subsubsection{Crank-Nicholson}
Ennél a módszernél az FTCS-et úgy módosítjuk, hogy a térbeli deriváltat az $n$-edik és az $n+1$-edik időpillanat
átlagából számoljuk, azaz:
\be
\frac{\rho_j^{n+1} - \rho_j^n}{\Delta t} +
v \frac{\rho_{j+1}^n - \rho_{j-1}^n + \rho_{j+1}^{n+1} - \rho_{j-1}^{n+1}}{4 \Delta x} = 0
\ee
Ennél a módszer amplitúdóban pontos, és feltétel nélkül stabil.
\subsection{Diffúziós egyenlet}
A diffúziós egyenlet a következőképpen néz ki:
\be
\frac{\partial T}{\partial t} = \alpha \frac{\partial^2 T}{\partial x^2}
\ee
Itt az
\be
s = \alpha \frac{\Delta t}{\Delta x^2}
\ee
változót vezetjük be.
\subsubsection{FTCS}
Ugye ez a módszer a Forward in Time, Centered in Space, és a diszkretizációja a következőképpen néz ki:
\be
\frac{T_j^{n+1}-T_j^n}{\Delta t} = \alpha \frac{T_{j+1}^n - 2T_j^{n} + T_{j-1}^n}{\Delta x^2}.
\ee
Ezt átrendezve:
\be
T_j^{n+1} = T_j^n + s (T_{j+1}^n - 2T_j^{n} + T_{j-1}^n)
\ee
Felírjuk a stabilitási egyenletet:
\be
G = 1 + s \left(\e^{\ii k \Delta x} - 2 + \e^{- \ii k \Delta x}\right)
\ee
Ezt átalakítva:
\be
G = 1 - 4 s \sin^2 \frac{k \Delta x}{2} < 1,
\ee
azaz ez a módszer az $s<1/2$ feltétellel stabil.
\subsubsection{Richardson}
Ez a módszer a kontinuitási egyenletnél említett Leapfrog módszernek (lásd \ref{sss:lf}) felel meg.
Itt az időbeli deriváltat is szimmetrikusan vesszük:
\be
\frac{T_j^{n+1}-T_j^{n-1}}{2 \Delta t} = \alpha \frac{T_{j+1}^n - 2T_j^{n} + T_{j-1}^n}{\Delta x^2}.
\ee
Ez a módszer itt azonban instabil.
\subsubsection{Du Fort -- Frankel}
Ez az előzőn javít, a térbeli deriváltat módosítjuk csak kicsit:
\be
\frac{T_j^{n+1}-T_j^{n-1}}{2 \Delta t} = \alpha \frac{T_{j+1}^n - (T_j^{n+1}+T_j^{n-1}) + T_{j-1}^n}{\Delta x^2}.
\ee
Ez a módszer feltétel nélkül stabil!
\subsubsection{Implicit módszer}
Ez ugyanaz, mint az FTCS, csak a térbeli deriváltat az ``új" időpillanatban vesszük:
\be
\frac{T_j^{n+1}-T_j^n}{\Delta t} = \alpha \frac{T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}}{\Delta x^2}.
\ee
Ez a módszer is feltétel nélkül stabil, ugyanakkor implicit, tehát egyenletrendszert kell megoldani az időben való
előrelépéshez. Ez megtehető különféle algoritmusokkal, pl a Thomas-algoritmus ilyen. Ezt most nem részletezném inkább.
\subsubsection{Crank-Nicholson}
Itt a térbeli deriváltat a két időpillanat átlagaként vesszük:
\be
\frac{T_j^{n+1}-T_j^n}{\Delta t} = \alpha \left[
\frac{T_{j+1}^{n+1} - 2T_j^{n+1} + T_{j-1}^{n+1}}{\Delta x^2} + \frac{T_{j+1}^n - 2T_j^{n} + T_{j-1}^n}{\Delta x^2}
\right].
\ee
Ez a módszer is feltétel nélkül stabil, és másodrendben pontos is.
\subsection{Transzport egyenlet}
A transzport egyenlet a kontinuitási és a diffúziós egyenletek keveréke:
\be
\frac{\partial T}{\partial t} + v \frac{\partial T}{\partial x} =
\alpha \frac{\partial^2 T}{\partial x^2}
\ee
Íme egy összefoglalás a fent említett módszerek ezen az egyenleten való alkalmazhatóságáról, azaz arról, hogy melyik
módszer mikor stabil és mikor pontos. Itt használjuk a korábban bevezetett $C$ és $s$ konstansokat.
\bigskip

\begin{center}
\begin{tabular}{||c||c||c||}
  \hhline{|t:=:t:=:t:=:t|}
  % after \\: \hline or \cline{col1-col2} \cline{col3-col4} ...
  Név & Stabilitás & Pontosság\\
  \hhline{|:=::=::=:|}
  FTCS & $C^2 < 2s < 1$ & $C^2 \ll 2s$\\ \hhline{||-||-||-|}
  Richardson-Leapfrog &$s=0$&  \\ \hhline{||-||-||-|}
  Du Fort -- Frankel & $|C| < 1$ &  $C^2 \ll 1$ \\ \hhline{||-||-||-|}
  Upwind & $|C|+2s<1$ & $\Delta t \Delta x \ll 2 / (1-|C|)$ \\ \hhline{||-||-||-|}
  Lax-Wendorff & $C^2+2s<1$ & $\Delta t \Delta x \ll 2$\\ \hhline{||-||-||-|}
  Crank-Nicholson & Stabil!! & $\Delta t \Delta x \ll 2$\\
  \hhline{|b:=:b:=:b:=:b|}
\end{tabular}
\end{center}
\section{Függelék}
\subsection{Paraméterek beolvasása fájlból}
Ha van egy\\
\verbß  egy 1ß\\
\verbß  ketto 2ß\\
\verbß  harom 3ß\\
tartalmú \verb"in.dat" fájlunk, és lefuttatjuk a\\
\verbß  #include <iostream>ß\\
\verbß  #include <string>ß\\
\verbß  #include <fstream>ß\\
\verbß  ß\\
\verbß  int main()ß\\
\verbß  {ß\\
\verbß    ifstream infile;ß\\
\verbß    infile.open("in.dat",ios::in);ß\\
\verbß    string nev;ß\\
\verbß    double ertek;ß\\
\verbß    while(infile >> nev >> ertek)ß\\
\verbß    {ß\\
\verbß      if(nev=="harom")ß\\
\verbß        cout << "Es a harmadik:" << endl;ß\\
\verbß      cout << nev << "=" << ertek << endl;ß\\
\verbß    }ß\\
\verbß    cout << "Ennyi." << endl;ß\\
\verbß    return 0;ß\\
\verbß  }ß\\
\verb"proba.cc" programot, akkor az a következőt írja ki:\\
\verbß  csanad@login03:~/temp> g++ proba.cc -o proba.exeß\\
\verbß  csanad@login03:~/temp> ./proba.exeß\\
\verbß  egy=1ß\\
\verbß  ketto=2ß\\
\verbß  Es a harmadik:ß\\
\verbß  harom=3ß\\
\verbß  Ennyi.ß\\
\verbß  csanad@login03:~/temp>ß\\
\subsection{Példa osztályok használatára}
Egy osztályt a következőképpen definiálhatunk:\\
\verbß  class diffegyß\\
\verbß  {ß\\
\verbß    private:ß\\
\verbß      int N;ß\\
\verbß      double a, b, x0, f0;ß\\
\verbßß\\
\verbß    public:ß\\
\verbß      string input, output;ß\\
\verbß      ifstream infile;ß\\
\verbß      ofstream outfile;ß\\
\verbß      diffegy(string input2, string output2);ß\\
\verbß      double func(double xx, double ff);ß\\
\verbß      double step(double xx, double ff, double ddx);ß\\
\verbß      void advance();ß\\
\verbß      ~diffegy();ß\\
\verbß  };ß\\
Ekkor lesznek belső (\verbßprivateß) és külső (\verbßpublicß) változói. Az osztály tetszőleges tagfüggvénye
(pl. a \verbßdouble func(double xx, double ff)ß) elérheti ezeket a változókat, azonban a belsőket \emph{csak}
a tagfüggvények érhetik el.

Az osztállyal azonos nevű függvény az osztály \emph{konstruktora}, az osztály egy objektumát inicializálva ez a függvény
fut le. A \verbß~diffegy()ß függvény az osztály destruktora, a program végén fut le.

Most definiáljuk a konstruktort! A két változója a ki- és bemeneti fájl neve, benne iniciálizáljuk a változók értékeit
és a két fájlt.\\
\verbß  diffegy::diffegy(string input2, string output2)ß\\
\verbß  {ß\\
\verbß    output=output2;ß\\
\verbß    input=input2;ß\\
\verbß    double ertek;ß\\
\verbß    string tipus;ß\\
\verbßß\\
\verbß    infile.open(input,ios::in);ß\\
\verbß    if(!infile)ß\\
\verbß    {ß\\
\verbß      cout<<input<<" nem megnyithato!"<<endl;ß\\
\verbß      exit(0);ß\\
\verbß    }ß\\
\verbß    elseß\\
\verbß      cout<<"Parameterek beolvasasa "<<input<<"-bol..."<<endl;ß\\
\verbß    while(infile>>tipus>>ertek)ß\\
\verbß    {ß\\
\verbß      if(tipus=="#a") a = ertek;ß\\
\verbß      if(tipus=="#b") b = ertek;ß\\
\verbß      if(tipus=="#x0") x0 = ertek;ß\\
\verbß      if(tipus=="#f0") f0 = ertek;ß\\
\verbß      if(tipus=="#N") N = (int)ertek;ß\\
\verbß      if(tipus=="#dx0") dx = ertek;ß\\
\verbß      if(tipus=="#df0") df = ertek;ß\\
\verbß    };ß\\
\verbß    infile.close();ß\\
\verbß    cout<<"Parameterek beolvasva."<<endl;ß\\
\verbßß\\
\verbß    outfile.open(output,ios::out);ß\\
\verbß    if(!outfile)ß\\
\verbß    {ß\\
\verbß      cout<<output<<"nem megnyithato!"<<endl;ß\\
\verbß      exit(0);ß\\
\verbß    }ß\\
\verbß    elseß\\
\verbß      cout<<output<<" megnyitasa irasra..."<<endl;ß\\
\verbßß\\
\verbß    outfile<<"#Kimeneti fajl. Oszlopok: x, f"<<endl;ß\\
\verbß    outfile.close();ß\\
\verbß  };ß\\
Ezek után definiáljuk a ``lépő" függvényt és a diffegyenlet jobb oldalát:\\
\verbß  double step(double xx, double ff, double ddx)ß\\
\verbß  {ß\\
\verbß    //ß pl Runge-Kutta lepes kerulhet ide, lásd \ref{p:rk4}\\
\verbß  };ß\\
\verbßß\\
\verbß  double func(double xx, double ff)ß\\
\verbß  {ß\\
\verbß    //ß itt kell definiálni a 'jobb oldalt', lásd \ref{p:adapt}\\
\verbß  };ß\\
Most pedig definiáljunk egy szubrutint, ami már a kimeneti fájlba írja az eredményt:\\
\verbß  void calculate()ß\\
\verbß  {ß\\
\verbß    /ß/ide kerul a korábbi fő-függvény, lásd \ref{p:adapt}\\
\verbß  };ß\\
Végül a destruktor:\\
\verbß  diffegy::~diffegy()ß\\
\verbß  {ß\\
\verbß    cout<<"Desctructing done."<<endl;ß\\
\verbß  };ß\\
És már jöhet is a fő-függvény. Ennek most van két argumentuma, azaz\\
\verbß  csanad@itl84:~/temp>./diffegy.exe input.dat output.datß\\
módon lehet meghívni. Ebben inicializálunk egy \verbßAß nevű, \verbßdiffegyß
típusú objektumot, majd meghívjuk ennek a \verbßcalculate()ß tagfüggvényét,
és kész is!\\
\verbß  int main(int argc, char *argv[])ß\\
\verbß  {ß\\
\verbß    int tip;ß\\
\verbß    if(argc!=3)ß\\
\verbß    {ß\\
\verbß      cerr<<"args: inputfile outputfile"<<endl;ß\\
\verbß      return 0;ß\\
\verbß    };ß\\
\verbßß\\
\verbß    cout<<"Inputfile: " <<argv[1]<<"; ";ß\\
\verbß    cout<<"Outputfile: "<<argv[2]<<endl;ß\\
\verbßß\\
\verbß    diffegy A(argv[1],argv[2]);ß\\
\verbß    cout<<"Constructing done."<<endl;ß\\
\verbßß\\
\verbß    A.calculate();ß\\
\verbßß\\
\verbß    return 0;ß\\
\verbß  };ß\\
\subsection{Animációk készítése gnuplottal}
Animációt a következőképpen készíthetünk gnuplot segítségével. Először is létrehozzuk az adatfájlt úgy, hogy az egyes
időpillanatoknak megfelelő adatok külön oszlopokba vannak írva. Ezek után készítünk egy plot fájlt, és elmentjük mondjuk
\verbßanimate.pltß néven:\\
\verbß  i=i+1ß\\
\verbß  plot "$0" u i w lß\\
\verbß  pause $1ß\\
\verbß  if(i<N) rereadß\\
Itt az \verbßiß változó fog futni \verbßNß-ig. A \verbß$0ß és \verbß$1ß változók az \verbßanimate.pltß program
betöltésekor adott argumentumok, a második azt határozza meg, hogy mekkora szünet legyen az egyes frame-ek között.
A fenti plot fájlt egy \verbßmain.pltß script fogja meghívni, ez a következőképpen néz ki:\\
\verbß  i=0ß\\
\verbß  N=80ß\\
\verbß  unset keyß\\
\verbß  set xrange [0:100]ß\\
\verbß  set yrange [0:20]ß\\
\verbß  call 'animate.plt' 'data.txt'ß\\
A fenti kis scriptben szereplő számokat természetesen szintén argumentumként is megadhatjuk, a \verbß$0ß, \verbß$1ß
\dots használatával. A scriptet gnuplotban a\\
\verbß  gnuplot> load 'main.plt'ß\\
paranccsal hívhatjuk meg.
Ha viszont a \verbßmain.pltß legelejére beírjuk a\\
\verbß  #!/usr/bin/gnuplotß\\
sort, akkor közvetlenül shellből is futtatható lesz (amennyiben a gnuplot valóban a fenti 
helyen van -- ez a \verbßwhich gnuplotß paranccsal ellenőrizhető).
\subsection{Oszlopok írása fájlba}
A \verbßC++ß nyelvben nincs közvetlen lehetőség arra, hogy egy fájl n.\ sorához hozzáfűzzünk
egy karaktert. Azonban létezik például egy \verbßgetline()ß parancs, amellyel ki lehet
olvasni a teljes aktuális sort egy karakter-állományba. Ehhez aztán hozzáfűzhetünk egy
számot, egy stringet, bármit, és beírhatjuk egy új fájlba.

Erre példa a következő szubrutin. Négy argumentuma van, az eredeti és az új fájl neve, 
illetve egy vektor, amelynek elemit a sorok végéhez hozzáfűzi, végül a sorok maximális
mérete, ez utóbbit tárolja a \verbßsizeß integer.
A szubrutin tehát:\\
\verbß  int append(char * inf, char * outf, double * v, int size)ß\\
\verbß  {ß\\
\verbß    char * temp;ß\\
\verbß    ifstream infile;         ß//a bemeneti fájl\\
\verbß    infile.open(inf,ios::in);ß\\
\verbß    ofstream outfile;        ß//a kimeneti fájl\\
\verbß    outfile.open(outf,ios::out);ß\\
\verbßß\\
\verbß    for(int i = 0; i < 5; i++)ß\\
\verbß    {ß\\
\verbß      temp = new char[size]; ß//ideiglenes karaktertömb\\
\verbß      infile.getline(temp,size);ß\\
\verbß      outfile << temp << setw(10) << v[i] << endl;ß\\
\verbß      delete temp;ß\\
\verbß    }ß\\
\verbß    infile.close();ß\\
\verbß    outfile.close();ß\\
\verbß    return 1;ß\\
\verbß  }ß\\
Ezt felhasználva a \verbßmain()ß függvénybe írhatjuk például a következőt:\\
\verbß  int main()ß\\
\verbß  {ß\\
\verbß    ofstream outfile;ß\\
\verbß    outfile.open("data.txt",ios::out);ß\\
\verbß    double * v;ß\\
\verbß    v = new double[5];ß\\
\verbß    for(int j = 0; j < 5; j++)ß\\
\verbß    {ß\\
\verbß      v[j] = 2.0*j + 2.0;ß\\
\verbß      outfile << setw(10) << v[j] << endl;ß\\
\verbß    }                        ß//itt beírtuk az első oszlopot,\\
\verbß                             ß//amely a kezdeti feltételt jelenti\\
\verbß    outfile.close();ß\\
\verbß    for(int j = 1; j < 6; j++)ß\\
\verbß    {ß\\
\verbß      calc(v);               ß//kiszámoljuk az új \verbßvß-t\\
\verbß      append("data.txt","temp.txt",v,2*j*11);ß\\
\verbß      calc(v);               ß//megint újra kiszámoljuk\\
\verbß      append("temp.txt","data.txt",v,(2*j+1)*11);ß\\
\verbß    }ß\\
\verbß    return 1;ß\\
\verbß  }ß\\

Ezt a programot lefuttatva, ha valahol definiáltuk a \verbßcalc(double * v)ß függvényt, akkor 13 oszlopot írunk a
\verbßdata.txtß fájlba egymás után, és csinálunk egy ideiglenes \verbßtemp.txtß fájlt is, amelyb az első 12 oszlopot
tartalmazza.

Természetesen ez csak egy nulladik lépés, egy példa, hogy hogyan lehet ilyesmit csinálni, hogyan lehet oszlopokat írni
egy fájlba anélkül, hogy a fájl teljes tartalmát előre ismernénk. A feladat konkrét jellegének megfelelően kell a fenti
kezdeményt módosítgatni, pl. dinamikusan lefoglalni a memóriában a tömböket, osztályokat használni. A fenti leírás
csak egy keretet, egy ötletet akar adni ehhez.
\end{document}
