A probléma: $L$ hosszúságú homogén fémrúd egyik végét $\tau_1$ hőmérsékleten tartjuk, másik végét $q$ hőárammal fűtjük. Kezdetben a rúd hőmérséklet-eloszlása lineáris, egyik végén a hőmérséklet $\tau_0$, a másik végén $\tau_1$. Határozzuk meg a hőmérséklet időbeli változását a rúd mentén!
A $T(x,t)$ hőmérsékletet keressük az $x$ koordináta és a $t$ idő függvényében.
Hővezetés esetén a hőáram-sűrűség (a felületegységen időegység alatt átáramlott hő) egy dimenzióban
$$j=-\kappa \frac{\partial T}{\partial x}\;.$$
Itt $\kappa$ a hővezetési együttható.
A rúd egy $x$ és $x+dx$ közötti szakaszának $dT$ hőmérsékletváltozásához
$$dQ=c(\rho A dx)dT$$
hőmennyiség közlése szükséges, ahol $c$ a fajhő (tömegegységre vonatkoztatva), $\rho$ a sűrűség és $A$ a keresztmetszet. A zárójelben levő szorzat a szóban forgó szakasz tömege. Másfelől a hőközlés a hővezetésből származik, a belépő és kimenő hőáramok különbségének megfelelően:
$$dQ=\left(-j(x+dx)+j(x)\right)A dt=\left(\left.\frac{\partial T}{\partial x}\right|_{x+dx}-\left.\frac{\partial T}{\partial x}\right|_{x}\right)\kappa A dt\approx \kappa A dt dx \frac{\partial^2 T}{\partial x^2}$$
A két kifejezést összevetve kapjuk a hővezetés egyenletét:
$$\frac{\partial T}{\partial t}=\frac{\kappa}{c \rho}\frac{\partial^2 T}{\partial x^2}\;.$$
Ezt a parciális differenciálegyenletet kell tehát megoldanunk a
$$T(L,t)=\tau_1$$
Dirichlet-féle határfeltétel, a
$$j(0,t)=-\kappa\left.\frac{\partial T}{\partial x}\right|_{0}=\frac{q}{A}$$
Neumann-féle határfeltétel, és a
$$T(x,0)=\tau_0\left(1-\frac{x}{L}\right)+\tau_1\frac{x}{L}$$
kezdeti feltétel esetében. Ezekkel a feltételekkel a feladat korrekt kitűzésű.
Az első lépés az egyenlet térbeli diszkretizálása. A problémát egydimenziósnak tekintjük, vagyis a keresztirányú méretek feltevés szerint jóval kisebbek a hosszirányúaknál. Felveszünk $N+1$ darab osztópontot $h=L/N$ távolságonként (tehát ekvidisztáns módon). Az egyes pontok $x$-koordinátái
$$x_j=jh\;,\quad j=0,1,2,\dots N$$
A diffúziós egyenletet a cellákra integráljuk:
$$\frac{d}{dt}\int_{x_j-h/2}^{x_j+h/2}T(x,t)dx=\frac{\kappa}{c
\rho}\left(\left.\frac{\partial T}{\partial
x}\right|_{x_j+h/2}-\left.\frac{\partial T}{\partial
x}\right|_{x_j-h/2}\right)$$
Ebből közelítőleg:
$$\frac{dT_j}{dt}=\frac{\kappa}{c
\rho h^2}\left(T_{j+1}-2T_j+T_{j-1}\right)\;,\quad j=1,2,\dots N-1$$
Itt $T_j=T(x_j)$. Közönséges differenciálegyenlet-rendszerhez jutottunk. A közelítés másodrendig pontos:
Helyettesítsük be $T(x,t)$ $x_j$ körüli Taylor-sorát, azaz a
$$T(x,t)=T_j+ \left.\frac{\partial T}{\partial
x}\right|_{x_j}(x-x_j)+\frac{1}{2}\left.\frac{\partial^2 T}{\partial
x^2}\right|_{x_j}(x-x_j)^2+\frac{1}{6}\left.\frac{\partial^3 T}{\partial
x^3}\right|_{x_j}(x-x_j)^3+\frac{1}{24}\left.\frac{\partial^4 T}{\partial
x^4}\right|_{x_j}(x-x_j)^4+\dots $$
kifejtést:
$$\frac{1}{h}\int_{x_j-h/2}^{x_j+h/2}T(x,t)dx=T_j+\frac{h^2}{8}\left.\frac{\partial^2 T}{\partial
x^2}\right|_{x_j}+\dots $$
$$\frac{T_{j+1}-2T_j+T_{j-1}}{h^2}=\left.\frac{\partial^2 T}{\partial
x^2}\right|_{x_j}+\frac{h^2}{12}\left.\frac{\partial^4 T}{\partial
x^4}\right|_{x_j}+\dots $$
$$\frac{1}{h}\left(\left.\frac{\partial T}{\partial
x}\right|_{x_j+h/2}-\left.\frac{\partial T}{\partial
x}\right|_{x_j-h/2}\right)=\left.\frac{\partial^2 T}{\partial
x^2}\right|_{x_j}+\frac{h^2}{24}\left.\frac{\partial^4 T}{\partial
x^4}\right|_{x_j}+\dots $$
A közelítés hibája tehát $h^2$ rendű.
A határfeltételek implementációja:
$x=L$ esetén Dirichlet-határfeltétel,
$$T_N=\tau_1\;.$$
$x=0$ esetén Neumann-határfeltétel,
$$\left.\frac{\partial T}{\partial
x}\right|_0=-\frac{q}{\kappa A}$$
Az utóbbi derivált kiszámítása másodrendű pontossággal:
$$T_1=T_0+\left.\frac{\partial T}{\partial
x}\right|_0h+\frac{1}{2}\left.\frac{\partial^2 T}{\partial
x^2}\right|_0h^2+\dots$$
$$T_2=T_0+\left.\frac{\partial T}{\partial
x}\right|_02h+\frac{1}{2}\left.\frac{\partial^2 T}{\partial
x^2}\right|_04h^2+\dots$$
Ebből
$$\left.\frac{\partial T}{\partial
x}\right|_0=\frac{4T_1-3T_0-T_2}{2h}+{\mathcal O}(h^2)$$
Az $x=0$-beli határfeltételből ezzel
$$T_0=\frac{4}{3}T_1-\frac{1}{3}T_2+\frac{2}{3}\frac{qh}{\kappa A}$$
Az időderivált implementálása C (és C++) nyelven:
Itt
$$Q=\frac{q}{\kappa A}$$
és
$$D=\frac{\kappa}{c
\rho }\;.$$
A következő feladat a közönséges differenciálegyenlet-rendszer megoldása. A
legegyszerűbb explicit Euler-módszerrel:
$$T_j^{(n+1)}=T_j^{(n)}+dt \left.\frac{\partial T^{(n)}}{\partial
t}\right|_{x_j}$$
azaz
Ha $dt$ elég nagy, akkor instabilitás lép fel, a megoldás abszolút
értékben végtelenhez tart. Vizsgáljuk meg ugyanis, hogy mi lesz a sorsa egy
(numerikus eredetű) kis perturbációnak!
Legyen
$$\frac{\partial{\bf \delta T}^{(n)}}{\partial
t}={\bf A}{\bf \delta T}^{(n)}$$
az időderivált linearizált változata (homogén lineáris, a konstans összeadandó
tagok kiesnek!).
Ekkor
$${\bf \delta T}^{(n+1)}=\left({\bf 1}+{\bf A}\right){\bf \delta T}^{(n)}=\left({\bf 1}+{\bf A}\right)^{n}{\bf \delta T}^{(1)}\;.$$
Legyen az ${\bf 1}+{\bf A}$ mátrix $k$-adik sajátértéke $\lambda_k$, a
hozzátartozó sajátvektora pedig $ {\bf s}^{(k)}$, azaz
$$ \left({\bf 1}+{\bf A}\right){\bf s}^{(k)}=\lambda_k{\bf s}^{(k)}\;.$$
A kezdeti perturbáció ${\bf \delta T}^{(1)}$ vektorát kifejtjük az $ {\bf s}^{(k)}$ sajátvektorok szerint:
$${\bf \delta T}^{(1)}=\sum_{k}c_k{\bf s}^{(k)}\;,$$
ezzel
$${\bf \delta T}^{(n+1)}=\sum_{k}c_k\lambda_k^n{\bf s}^{(k)}\;.$$
Ha tehát
$$\left|\lambda_k\right|>1$$
valamelyik $k$ esetén, akkor a hiba az időben exponenciálisan növekedve végtelenhez tart. Jelen esetben
az ${\bf A}$ mátrix diagonális és mellékdiagonális elemei különböznek csak a
nullától (tridiagonális mátrix):
$$A_{jj}=-2\frac{D}{h^2} dt\;,\quad A_{jj\pm 1}=\frac{D}{h^2} dt$$
A sajátvektorokat
$$s_j^{(k)}=\exp\left(ipj\right)$$
képzetes kitevőjű exponenciálisként keressük. Ekkor
${\bf 1}+{\bf A}$ megfelelő sajátértékére
$$\lambda=1-4\frac{Ddt}{h^2}\sin^2p$$
adódik. (A $p$ hullámszám valójában diszkrét, és lehetséges értékei a
peremfeltételek alapján határozhatók meg, de ezzel most
nem foglalkozunk.)
Instabilitás akkor lép fel, ha
$$dt>\frac{h^2}{2Ddt}\;,$$
ekkor ugyanis $\lambda<-1$ ($p=\pi/2$ esetén).
