W celu przeprowadzenia symulacji należy określić algorytm rozwiązania matematycznego opisanego wyżej równania Fourier'a.

Paraboliczne równania różniczkowe cząstkowe

Zagadnienia transportu ciepła należą to klasy problemów, których analiza wymaga rozwiązania złożonych równań różniczkowych cząstkowych. Tylko w niektórych, bardzo uproszczonych przypadkach istnieją ścisłe rozwiązania analityczne, dlatego też do badania problemów przepływu ciepła, dyfuzji, czy też emisji i absorpcji promieniowania cieplnego coraz częściej stosowane są komputerowe metody numeryczne.  Chociaż ich właściwą wagę i znaczenie można docenić dopiero dla skomplikowanych układów, to jednak poznanie tych metod dla przykładów stosunkowo prostych, ma duże znaczenie dydaktyczne.

Często spotykane równanie różniczkowe cząstkowe występujące w wielu sytuacjach fizycznych ma postać:

gdzie D to stała dyfuzji (niekiedy zależna od położenia), a S to funkcja źródeł. Najbardziej znanym przykładem takiego równania jest równanie Schroedingera.

Problem numerycznego rozwiązania parabolicznych równań cząstkowych często sprowadza się do tzw. problemu wartości początkowej: mamy pole f w czasie początkowym, szukamy w czasie późniejszym, a ewolucja systemu podlega pewnym warunkom brzegowym (i tak np. funkcja gęstości prawdopodobieństwa zanika w dużej odległości, a temperatura lub przepływ ciepła jest określony na jakiejś powierzchni. Metody rozwiązywania takich problemów na komputerze są generalnie proste koncepcyjnie, chociaż trzeba mieć świadomo ść kilku subtelności (np. słaba stabilność metody prostej dyskretyzacji).

Dla przypadku jednowymiarowego przepływu ciepła w ciałach stałych równanie Fouriera przyjmuje postać:

zakładamy warunki brzegowe Dirichleta, czyli określamy wartoś ci pola na końcach rozpatrywanego przedziału - grzejnik i chłodnica.
Standardowe podejście analizy numerycznej do podobnych problemów przewiduje aproksymację różniczek przestrzennych skończnymi elementami różnicowymi na jednorodnej siatce składającej się z N+1 punktów o odstępie h = 1/N oraz różniczki cząstkowej formułą różniczkową pierwszego rzędu o kroku czasowym Dt. [13, 6]

 Metody różnic skończonych a metody relaksacyjne

Metoda różnic skończonych należy do najczęściej stosowanych metod numerycznych rozwiązywania równań różniczkowych cząstkowych. Polega ona na przechodzeniu w równaniu różniczkowym od pochodnych do odpowiednich ilorazów różnicowych - czyli przejściu od równań różniczkowych do równań różnicowych, które wiążą ze sobą wartości szukanej funkcji w pojedynczych, odosobnionych punktach. Punkty te są natomiast wybierane tak, aby tworzyły siatkę regularną (kwadratową, prostokątną lub sześcienną), a jej rodzaj zależy zwykle od rodzaju układu współrzędnych naturalnego dla badanego zagadnienia.

Wprowadzając zależność wartości szukanej funkcji od punktów sąsiednich w kroku czasowym t-1 można doprowadzić problem do rozwiązania n równań liniowych z wieloma niewiadomymi. Przy dużej liczbie węzłów n rozwiązanie otrzymanego układu n równań nawet za pomocą krakowianów, może być bardzo uciążliwe. Dlatego w praktyce stosuje się często tzw. metodę relaksacyjną. Polega ona na:

• sporządzeniu siatki kwadratowej z węzłami odległymi o h0 i przypisaniu wartoś ci węzłom brzegowym zgodnie z warunkiem brzegowym, natomiast węzłom wewnętrznym przypisujemy dowolne liczby;
• węzłom wewnętrznym przypisujemy następnie wartości obliczone ze wzoru różnicowego (uwzględniając węzły sąsiednie);
• czynność tę powtarzamy wielokrotnie (im więcej razy, tym większą dokładność otrzymamy)

Ciągi wartości szukanej funkcji w punktach węzłowych otrzymane metodą relaksacyjną są zbieżne do pierwiastków układu równań powstałego przy dyskretyzacji równania różniczkowego. Opis metod relaksacyjnych można znaleźć np. pracach [6, 9, 14].

Metoda prostej dyskretyzacji - stabilność metod różnic skończonych dla równań parabolicznych

Istnieje wiele metod różnic skończonych. Najistotniejszym czynnikiem różniącym te metody jest jawność zastosowanego schematu różnicowego.
Wybór jednej z metod różnicowych do rozwiązania danego problemu jest tylko pierwszym krokiem na drodze budowy całego algorytmu numerycznego. Metody ró żnic skończonych choć proste i dobrze określone charakteryzują się dużą ogólnością podejścia do zagadnień numerycznych. Dla całkowitego sprecyzowania stosowanego algorytmu trzeba jeszcze dokonać wyboru w jaki sposób różniczki są przybliżane różnicami skończonymi. Takich sposobów jest wiele, a dobry wybór najbardziej odpowiedniego z nich jest najtrudniejszą rzeczą przy numerycznym traktowaniu podobnyc h zagadnień. Od tego jakie podejście zostanie przyjęte zależy przede wszystkim stabilność metody numerycznej, dokładność, zbieżność, a co za tym idzie ogólna przydatność algorytmu. Bogaty opis różnych metody znajduje się w pracy [6].

Można zastosować najprostsze rozwiązanie i przybliżyć różniczki przez różnice skończone pierwszego rzędu. Zakładając, że mamy jednorodną siatkę N+1 punktów o odstępie h=1/N. Niech n oznacza indeks określający krok czasowy. Wtedy równanie Fouriera można przepisać w postaci:

pozwala to już na zdefiniowanie schematu różnicowego. Biorąc pod uwagę, że drugą pochodną przestrzenną można przedstawić w postaci:

schemat różnicowy przybiera formę:

(Hf)i = 1/h2 * (d2f)i

(3.6)

co w efekcie daje poszukiwaną formułę zależności funkcji f w puktach kroku czasowego n+1 od wartości w kroku czasowym n.

Tak określona formuła różniczkowa w sposób poprawny definiuje rozwój układu fizycznego. Wymaga ona jednak przyjęcia bardzo małego kroku czasowego (prawie dwa rzędy wielkości mniejszego niż charakterystyczny krok czasowy rozpatrywanego modelu Dt < h²/2) [13] Kryterium stabilno ści jawnego schematu różnicowego oraz dokładniejsze oszacowanie wielkości kroku czasowego można znaleźć w pracy [6].  Źródło dużej wrażliwości formuły na wielkość kroku czasowego leży w samej metodzie. Każdy profil rozwiązania równania można przedstawić w postaci szeregu Fouriera, każdy element tego szeregu spełnia równanie i jest jego rozwiązaniem szczególnym. Funkcje trygonometryczne można przedstawić w postaci exponencjalnej, a funkcję exp przybliżyć różnicowo z rożną dokładnością. Przyjmując pierwszy rząd rozwinięcia exp w szereg zawężamy znacznie możliwość operowania krokiem czasowym.

Metoda ulepszonej dyskretyzacji - metoda niejawna

Najprostszą metodą uniknięcia niestabilności pr ostej dyskretyzacji równania jest zastąpienie różniczki przestrzennej drugiego rzędu branej w czasie n (starym), przez tę w czasie n+1 (nowym).

Schemat różnicowy będzie wyglądał następująco:

rozwiązanie tego równania ze względu na fⁿ⁺¹prowadzi do równania operatorowego

dla bardzo małego kroku czasowego jest ono równoważne podejściu poprzedniemu, jednak dla większych Dt nie traci stabilności numerycznej.

Można nawet zastosować dokładniejszą formułę - przyjmując za różniczkę przestrzenną średnią ze schematów różnicowych dla dwóch sąsiednich kroków czasowych.

czyli

Z punktu widzenia wydajności metody ulepszonej dyskretyzacji ujemną jej cechą jest konieczność znalezienia w każdym kroku czasowym rozwiązania zbioru liniowych równań. Jest to równoznaczne podziałaniu odwróconą macierzą  (1 + HDt/2)  na wektor z nawiasów kwadratowych (3.11). Ponieważ elementy macierzy nie zależą od poszczególnych kroków czasowych można oszczędzić na obliczeniach i wyznaczyć macierz odwrotną jeden tylko raz na samym początku algorytmu. Aplikacja nadal jeszcze wymaga N² operacji jeśli obliczenia mają być przeprowadzane wprost.

Metoda eliminacji Gaussa ze zwrotnym podstawianiem

Tylko w przypadku metody prostej dyskretyzacji mamy do czynienia z jawną formułą zależności wartości szukanej funkcji w kroku n+1 od jej wartości w kroku n. Ulepszenie algorytmu poprzez użycie jakiejkolwiek kombinacji wartości funkcji w innych krokach prowadzi z jednej strony do większej stabilności numerycznej algorytmu, ale z drugiej utrudnia znalezienie wartości szukanej funkcji. Formuła jest w takim wypadku niejawna, a rozwiązanie takiego problemu sprowadza się do odwrócenia macierzy trójprzekątniowej (czyli rozwiązania układu równań na wartości funkcji w różnych krokach przestrzennych).

Jedną z najskuteczniejszych metod znalezienia macierzy odwrotnej jest metoda eliminacji Gaussa, która lekko zmodyfikowana (uproszczona) dla macierzy trójprzekątniowej prowadzi do algorytmu rzędu N.

Dla ustalenia uwagi skoncentrujmy się teraz nad równaniem różnicowym:

Zadanie znalezienia nowej wartości f sprowadza się do obliczenia wartości własnych operatora 1+HDt.

Niech A_i = 1 + HDt oraz B_i = fⁿ+ S(n)Dt
W naszym pr oblemie otrzymujemy macierz trójdiagonalną :

W celu rozwiązania tak powstałego układu równań zakładamy, że szukane rozwiązanie spełnia rekurencyjnie zależność f_i+1=a_if_i + b_i,  a i b trzeba więc określić.  Podstawiając tę zależność do równania definiującego naszą macierz trójdiagonalną oraz rozwiązując f_i   otrzymujemy rekurencyjne zależności na a i b .

Należy więc wyznaczyć a i b zaczynając od końca siatki. Zwykle siatka jest numerowana od 0 do N-1, nasze obliczenia należy rozpocząć od węzła N-2 pamiętając o tym, żeby poprawnie uwzględnić warunki brzegowe.
a_N-1  = 0; b_N-1  = f_N

  Mając a i b wyznaczone dla całej siatki można obliczyć szukaną funkcję  w pętli rozpoczynającej się od 0.