Metody numeryczne fizyki/Rozwiązywanie równań różniczkowych cząstkowych z warunkami początkowymi

Rozwiązywanie równań różniczkowych cząstkowych z warunkami początkowymi

---

Licencja

Autor: Mirosław Makowiecki Absolwent UMCS Fizyki Komputerowej Uniwersytetu Marii Curie-Skłodowskiej w Lublinie Email: miroslaw(kropka)makowiecki(małpa)gmail(kropka)pl Dotyczy: książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami. Użytkownika książki, do której należy ta strona, oraz w niej zawartych stron i w nich podstron, a także w nich kolumn, wraz z zawartościami nie zwalnia z odpowiedzialności prawnoautorskiej nieprzeczytanie warunków licencjonowania. Umowa prawna: Creative Commons: uznanie autorstwa, na tych samych warunkach, z możliwością obowiązywania dodatkowych ograniczeń. Autor tej książki dołożył wszelką staranność, aby informacje zawarte w książce były poprawne i najwyższej jakości, jednakże nie udzielana jest żadna gwarancja, czy też rękojma. Autor nie jest odpowiedzialny za wykorzystanie informacji zawarte w książce, nawet jeśli wywołaby jakąś szkodę, straty w zyskach, zastoju w prowadzeniu firmy, przedsiębiorstwa lub spółki bądź utraty informacji, niezależnie czy autor (a nawet Wikibooks) został powiadomiony o możliwości wystąpienie szkód. Informacje zawarte w książce mogą być wykorzystane tylko na własną odpowiedzialność.

---

---

Będziemy się tutaj zajmowali równaniami cząstkowymi względem różnych parametrów wraz z warunkami brzegowymi.

Wprowadzenie do równań różniczkowych cząstkowych

Przykładem równania brzegowego jest:

$Lu={{\partial u} \over {\partial t}}-{{\partial ^{2}u} \over {\partial x^{2}}}=f(x,t)\;$

(8.1)

Rozwiązaniem równania różniczkowego (8.1) jest u=u(x,t), które jest spełnione wewnątrz prostokąta D={(x,t):0≤x≤l,0≤t≤T}, którego warunkami brzegowymi są u(x,0)=g(x), 0≤x≤l, a także u(0,t)=φ(t) oraz u(l,t)=ψ(t) dla 0≤t≤T. Opisane powyżej zagadnienia nazywamy zagadnieniami brzegowymi pierwszego rodzaju lub zagadnieniami Dirichleta.

Sprowadzanie pochodnych cząstkowych do zagadnień różnicowych

Weźmy sobie prostokąt o osiach położenia x i czasu t, który podzielimy na malutkie prostokąciki o kroku h względem osi położenia i o kroku k względem osi czasowej:

$x=ih{\mbox{, }}t=jk{\mbox{, }}i,j\in C\;$

(8.2)

Obierzmy sobie operator Λ, który w wyniku działania na funkcję u(x,t), otrzymamy pewne wyrażenie różnicowe, które możemy przedstawić względem pochodnych cząstkowych pierwszego rzędu względem czasu i drugiego rzędu względem położenia.

$\mathbf {\Lambda } u(x,t)={{u(x+h,t)-2u(x,t)+u(x-h,t)} \over {h^{2}}}=\;$
$={{\partial ^{2}u(x,t)} \over {\partial x^{2}}}-{{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+{{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+O(h^{2}+k^{2})=\;$
$={{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}-\left({{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}-{{\partial ^{2}u(x,t)} \over {\partial x^{2}}}\right)+O(h^{2}+k^{2})=\;$
$={{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}-{{k} \over {2}}{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}+O(h^{2}+k^{2})\;$

(8.3)

$\mathbf {\Lambda } u(x,t+k)={{u(x+h,t+k)-2u(x,t+k)+u(x-h,t-k)} \over {h^{2}}}=\;$
$={{\partial ^{2}u(x,t+k)} \over {\partial x^{2}}}-{{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+{{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+O(h^{2}+k^{2})=\;$
$={{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+\left({{\partial ^{2}u(x,t+k)} \over {\partial x^{2}}}-{{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}\right)+O(h^{2}+k^{2})=\;$
$={{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+{{k} \over {2}}{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}+O(h^{2}+k^{2})\;$

(8.4)

Weźmy teraz rodzinę operatorów różnicowych $L_{h,k}^{\sigma }\;$, który jest zależny od parametru σ i operatora Λ i na samym końcu od funkcji u(x,t), i na podstawie (8.3) i (8.4), mamy:

$L_{h,k}^{\sigma }u(x,t)={{u(x,t+k)-u(x,t)} \over {k}}-\left(\sigma \mathbf {\Lambda } u(x,t+k)+(1-\sigma )\mathbf {\Lambda } u(x,t)\right)=\;$
$={{\partial u(x,t+k/2)} \over {\partial t}}-{\Bigg [}\sigma \left({{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}+{{k} \over {2}}{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}\right)-\;$
$-\left(1-\sigma \right)\left({{\partial ^{2}u(x,t+k/2)} \over {\partial x^{2}}}-{{k} \over {2}}{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}\right){\Bigg ]}=\;$
$={{\partial u(x,t+k/2)} \over {\partial t}}+\left({{1} \over {2}}-\sigma \right)k{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}+O(h^{2}+k^{2})=\;$
$=Lu(x,t+k/2)+\left({{1} \over {2}}-\sigma \right)k{{\partial ^{3}u(x,t+k/2)} \over {\partial x^{2}\partial t}}+O(h^{2}+k^{2})\;$

(8.5)

Oznaczmy teraz przez u(x,t) dokładne rozwiązanie pewnego zagadnienia różnicowego, a v(x,t) rozwiązanie jego odpowiednika różnicowego. Rozpocznijmy rozważania nad warunkami brzegowymi, które to możemy przedstawić jako: v_i,0=g_i dla i=1,...,N-1, gdzie g_i=g(ih). W wezłach (0,j) i (N,j) możemy przedstawić warunki brzegowe v_0,j=φ_j i v_N,j=ψ_j dla j=0,1,...,M. W węzłach (i,j) podczas działania operatora $L_{h,k}^{\sigma }\;$ na funkcję v_ij zapisaną w kolejnych węzłach otrzymujemy równanie zależne od parametru σ oraz od "h", "k", "i" i "j":

$L_{h,k}^{\sigma }v_{ij}=f_{i,j+1/2}{\mbox{ dla }}i=1,...,N-1{\mbox{ oraz }}j\geq 1\;$

(8.6)

We wzorze (8.6) użyto oznaczenia f_i,j+1/2=f(ih,jk+k/2). Wykorzystując definicję $L_{h,k}^{\sigma }\;$ zapisaną na początku obliczeń (8.5), to równanie (8.6), przy założeniu α=k/h², możemy przepisać do postaci:

$k^{-1}\left(u_{i,j+1}-u_{i,j}\right)-\;$
$-\left[\sigma h^{-2}\left(u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}\right)+(1-\sigma )h^{-2}\left(u_{i+1,j}-2u_{i,j}+u_{u-1,j}\right)\right]=\;$
$=f_{i,j+1/2}\Rightarrow -\alpha \sigma v_{i+1,j+1}+(1+2\alpha \sigma )v_{i,j+1}-\alpha \sigma v_{i-1,j+1}=\alpha (1-\sigma )v_{i+1,j}+\;$
$+\left[1-2\alpha (1-\sigma )\right]v_{i,j}+\alpha (1-\sigma )v_{i-1,j}+kf_{i,j+1/2}$

(8.7)

W równaniu (8.7) należy przyjąć σ≠0, wtedy ono wiąże wartości v_ij w sześciu węzłach, ale już dla σ=0 lub σ=1 wartości v_ij wiążą w czterech węzłach. Gdy σ=0 równanie (8.7) możemy przepisać jako:

$v_{i,j+1}=\alpha v_{i+1,j}+(1-2\alpha )v_{i,j}+\alpha v_{i-1,j}+f_{i,j+1/2}\;$

(8.8)

Natomiast, gdy α≠0 równanie (8.7) zapisujemy w postaci macierzowej:

$\mathbf {A} \mathbf {V} _{j+1}=\mathbf {F} _{j+1}\;$

(8.9)

Co rozpisując poszczególne macierze i wektory występujące w równaniu (8.9), otrzymujemy:

$\mathbf {A} ={\begin{bmatrix}(1+2\alpha \sigma )&-\alpha \sigma &\\-\alpha \sigma &(1+2\alpha \sigma )&-\alpha \sigma &&0&\\&\ddots &\ddots &\ddots \\&0&\ddots &\ddots &-\alpha \sigma \\&&&-\alpha \sigma &(1+2\alpha \sigma )\end{bmatrix}}\;$

(8.10)

$\mathbf {V} _{j+1}={\begin{bmatrix}v_{1,j+1}\\v_{2,j+1}\\\vdots \\v_{N-2,j+1}\\v_{N-1,j+1}\end{bmatrix}}\;$

(8.11)

Macierz (8.10) jest macierzą silnie diagonalnie dominującą dla σ>0, zatem układy równań o takich macierzach mają jednoznaczne rozwiązanie. Ale już dla σ<0, macierz (8.10) już nie jest macierzą diagonalnie dominującą. Można udowodnić, że układy o takich macierzach dla pewnej wartości σ<0 macierz (8.10) może być macierzą osobliwą Oznaczmy przez $C_{n}^{m}({\overline {\mathbf {D} }})\;$, która oznacza zbiór funkcji określonych na zbiorze ${\overline {\mathbf {D} }}\;$, które ma m ciągłych pochodnych względem x i n ciągłych pochodnych względem t. Patrząc na obliczenia (8.5) dochodzimy do dwóch wniosków jedno dla σ=1/2 (dla $u\in C_{3}^{4}({\overline {\mathbf {D} }})\;$), a drugie dla σ≠1/2 (dla $u\in C_{2}^{4}({\overline {\mathbf {D} }})\;$):

$L_{h,k}^{1/2}u_{ij}=f_{i,j+1/2}+O(h^{2}+k^{2})\;$

(8.12)

$L_{h,k}^{\sigma }u_{ij}=f_{i,j+1/2}+O(h^{2}+k)\;$

(8.13)

Stabilność rozwiązań równań cząstkowych z warunkami brzegowymi i jego zbieżność

Weźmy dwie takie funkcje f(x,t), dwie funkcje g(x), która jest warunkiem brzegowym początkowym dla t=0 dla funkcji v(x,t), i na samym końcu dwie funkcje v. Jedne funkcje oznaczymy bez daszka, a drugie z daszkiem, wtedy możemy napisać oznaczenia: $V=v-{\tilde {v}}\;$, $F=f-{\tilde {f}}\;$ i $G=g-{\tilde {g}}\;$, wtedy równanie (8.7) możemy przepisać jako:

$-\alpha \sigma V_{i+1,j+1}+(1+2\alpha \sigma )V_{i,j+1}-\alpha \sigma V_{i-1,j+1}=\;$
$=\alpha (1-\sigma )V_{i+1,j}+[1-2\alpha (1-\sigma )]V_{i,j}+\alpha (1-\sigma )V_{i-1,j}+kF_{i,j+1/2}\;$
$V_{i,0}=G_{i}{\mbox{ dla }}i=1,2,..,N-1\;$

$V_{0,j}=V_{N,j}=0{\mbox{ dla }}j=0,1,...,M\;$

$F_{i,j+1/2}=F(ih,jk+k/2){\mbox{, }}G_{i}=G(ih)\;$

(8.14)

Warunki początkowe na równanie (8.12) są bardzo podobne do warunków początkowych (8.6), i dlatego nie będziemy tutaj rozważać. Warunki brzegowe jakie wchodzą do równania (8.14) zwykle są obarczone pewnymi błędami. Zaburzenia spowodowane funkcjami f,g,φ,ψ, które występują w równaniu (8.1) i w jego warunkach brzegowych powodują, że przy zagęszczaniu siatki powodują, że zmiany w poszczególnych warstwach wzrastają nieograniczenie. Tutaj będziemy dlatego rozważać, że f i g są obarczone pewnymi błędami, a pozostałe funkcje są są wyznaczane dokładnie, tzn. φ i ψ. Jeśli chcemy określić jakie są zmiany w kolejnych warstwach wprowadza się normę ||⋅|| w klasie funkcji, którego argumentami jest zmienna x_i=ih. Weźmy sobie Y=Y(x,t), wtedy przez Y^r oznaczmy funkcję, dla której dziedziną jest zbiór zmiennych oznaczony przez x_i, która przyjmuje wartości Y(ih,rk) dla t=rk. Wtedy możemy określić, czy (8.12) jest stabilna względem prawej strony i warunku początkowego, wtedy musimy określić stałe dodatnie, jeśli istnieją takie, że dla dowolnych funkcji G możemy napisać warunek:

$||V^{j}||\leq A_{1}||G||+A_{2}\cdot \max _{0\leq s\leq s-1}||F^{s+1/2}||{\mbox{ dla }}j=1,...,M\;$

(8.15)

Jeśli dodatkowo mamy tożsamościowo F=0 lub G=0, to wtedy na podstawie tego dla zachodzącej tożsamości (8.14) możemy określić metodę, którą nazwiemy stabilną względem warunku początkowego lub stabilną względem prawej strony, która wedle (8.13) możemy zapisać w postaci dwóch nierówności :

$||V^{j}||\leq A_{1}||G||\;$

(8.16)

$||V^{j}||\leq A_{2}\cdot \max _{0\leq s\leq j-1}||F^{s+1/2}||{\mbox{ dla }}j=1,...,M\;$

(8.17)

Widzimy, że małe zmiany funkcji g(x) dla t=0, a także małe zmiany funkcji f dla t=(j+1/2)k, dla j=0,1,..,M-1 powodują małe zmiany rozwiązania całego funkcji u(x,t). Zdefiniujmy normę:

$||V||=||V||_{c}=\max _{0\leq i\leq N}|V(x_{i})|\;$

(8.18)

Naszym dalszym celem jest udowodnienie nierówności, którą napiszemy w postaci poniżej zależne od normy ||G|| oraz od ||F^j+1/2||:

$||V^{j}||_{c}\leq ||G||_{c}+t_{j}\cdot \max ||F^{i+1/2}||{\mbox{ dla }}j=1,...,M{\mbox{, gdy }}\sigma \geq 1-{{h^{2}} \over {2k}}{\mbox{ dla }}0\leq \sigma \leq 1\;$

(8.19)

Określmy, teraz dwie funkcje y i Y, które są określone na zbiorze x_i dla i=0,1,..,N, wtedy zdefiniujmy taki operator P, który na działaniu na funkcję y_i otrzymujemy pewną kombinację liniową y_i+1, y_i, y_i-1:

$Py_{i}=\alpha \sigma y_{i+1}-(1+2\alpha \sigma )y_{i}+\alpha \sigma y_{i-1}{\mbox{, }}i=1,...,N-1{\mbox{, }}\alpha \sigma >0\;$

(8.20)

Twierdzenie pierwsze
Weźmy sobie ciąg niestałej funkcji y i jeśli zachodzi Py_i≥0 (Py_i≤0), to wtedy y_i nie może mieć wartości największej dodatniej (najmniejszej ujemnej) w punktach x_i określonych dla i=1,..,N

Dowód
Udowodnijmy pierwszą część twierdzenia powyższego, gdy Py_i≥0, to y_i nie może przyjmować wartości największej dodatniej. Będziemy dowodzili to twierdzenie nie wprost. W tym celu weźmy takie największe y_s, by by przyjmowało wartość $y_{s}=\max _{0\leq i\leq N}y_{i}=M>0\;$, wtedy napiszmy prawą stronę (8.18) dla i=s w sposób:

$Py_{s}=\alpha \sigma y_{s+1}-(1+2\alpha \sigma )y_{s}+\alpha \sigma y_{s-1}<\alpha \sigma M-(1-2\alpha \sigma )M+\alpha \sigma M=-M<0\;$

(8.21)

Otrzymaliśmy w (8.21), że to wyrażenie Py_i jest mniejsze od zera, a powinno być większe lub równe niż zero, co otrzymaliśmy sprzeczność, co znaczy, że udowodniliśmy pierwszą część powyższego twierdzenia. Drugą jego część dowodzimy podobnie.

Twierdzenie drugie
Jeśli zachodzi nierówność pierwsza |Py_i|≤-PY_i dla i=1,...,N-1, nierówność druga |y₀|≤Y₀ i na ostatku nierówność trzecia |y_n|≤Y_N, to wtedy zachodzi |y_i|≥Y_i dla i=0,1,..,N.

Dowód
Z nierówności pierwszej dowodzonego twierdzenia możemy powiedzieć P(Y_i-y_i)<0 i P(Y_i+y_i)<0, co z twierdzenia pierwszego wynika, że wartości Y_i-y_i, Y_i+y_i przyjmują wartości dodatnie, stad biorąc jeszcze drugą i trzecią nierówność dowodzonego twierdzenie udowodniliśmy, że zachodzi |y_i|≤Y_i

Twierdzenie trzecie
Jeśli zachodzi PY_i=-H dla H≥0, a także Y₀=Y_N=0, to zachodzi nierówność Y_i≤H dla i=0,1,...,N

Dowód
Dla i=0 i i=N na pewno zachodzi Y_i≤H, By udowodnić nasze twierdzenie dla i=1,..,N-1 weźmy takie Y_s, dla s=1,..,N-1, by przyjmowało wartość największą równą H, wtedy możemy napisać:

$PY_{s}=\alpha \sigma Y_{s+1}-(1+2\alpha \sigma )Y_{s}+\alpha \sigma Y_{s-1}=-\alpha \sigma (Y_{s+1}-Y_{s})-\alpha \sigma (Y_{s}-Y_{s-1})-Y_{s}\leq -Y_{s}=-H\;$

(8.22)

Stąd dochodzimy do wniosku, że Y_i≤H dla i=0,1,..,N.

Wzór (8.12) możemy zapisać, wiedząc, że istnieje największa wartość V_i,j i F_i,j+1/2 po prawej stronie naszego równania:

$\alpha \sigma V_{i+1,j+1}-(1+2\alpha \sigma )V_{i,j+1}+\alpha \sigma V_{i-1,j+1}=-\max _{0\leq i\leq N}|V_{i,j}|-k\max _{0\leq iN}|F_{i,j+1/2}|=-||V^{j}||-k||F^{j+1/2}||\;$

(8.23)

Obierzmy sobie ograniczenie na σ, tzn. pierwsze ograniczenie 0≤σ≤1, i drugie ograniczenie wykorzystując definicję α, wtedy:

$\sigma \geq 1-{{h^{2}} \over {2k}}=1-{{1} \over {2\alpha }}\Rightarrow 2\sigma \alpha \geq 2\alpha -1\Rightarrow 1+2\alpha \sigma \geq 2\alpha \;$

(8.24)

Jesli wybierzemy największą wartość V_i,j względem "i" i wykorzystując (8.22), wtedy otrzymujemy (8.23) i pamiętając o ograniczeniu na σ, dalej używając jeszcze twierdzenie drugie i trzecie, dochodzimy do wniosku:

$||V^{j+1}||\leq ||V^{j}||+k||F^{j+1/2}||\;$

(8.25)

Pierwszy wyraz w (8.23) po jego prawej stronie możemy rozwinąć używając tego samego równania aż do ||V⁰||, otrzymujemy:

$||V^{j+1}||_{c}\leq ||V^{j}||_{c}+k||F^{j+1/2}||_{c}\leq ||V^{j-1}||_{c}+||F^{j-1+1/2}||_{c}+||F^{j+1/2}||_{c}\leq \cdots \;$
$\leq ||V^{0}||_{c}+k\sum _{s=0}^{j}||F^{s+1/2}||\leq ||G||_{c}+k(j+1)\max _{0\leq s\leq j}||F^{s+1/2}||\Rightarrow \;$
$\Rightarrow ||V^{j}||_{c}\leq ||G||_{c}+kj\max _{0\leq s\leq j-1}||F^{s+1/2}||=||G||_{c}+T\max _{0\leq s\leq j-1}||F^{s+1/2}||\;$

(8.26)

Weźmy sobie równanie (8.14) i do niego podstawmy rozwiązanie zależne od λ_s, która nie zależy od "i" i "j", które zapisujemy w postaci:

$V_{ij}^{(s)}=\lambda _{s}^{j}\sin {{s\pi i} \over {N}}{\mbox{, }}s=1,\cdots N-1\;$

(8.27)

Co po podstawieniu (8.27) otrzymujemy:

$-\alpha \sigma \lambda _{s}\left[\sin {{s\pi (i+1)} \over {N}}-2\sin {{s\pi i} \over {N}}+\sin {{s\pi (i-1)} \over {N}}\right]+\lambda _{s}\sin {{s\pi } \over {N}}=\;$
$=\alpha (1-\sigma )\left[\sin {{s\pi (i+1)} \over {N}}-2\sin {{s\pi i} \over {N}}+\sin {{s\pi (i-1)} \over {N}}\right]+\sin {{s\pi i} \over {N}}\;$

(8.28)

We wzorze (8.28) będziemy wykorzystywać wzory z trygonometrii, co tuż po wykorzystaniu nasze równanie dzielimy przez $_{\sin {{s\pi i} \over {N}}}\;$, w ten sposób otrzymujemy tożsamość:

$-\alpha \sigma \lambda _{s}\left[2\cos {{s\pi } \over {N}}-2\right]+\lambda _{s}=\alpha (1-\sigma )\left[2\cos {{s\pi } \over {N}}-2\right]+1\Rightarrow \lambda _{s}={{1-\alpha (1-\sigma )4\sin ^{2}{{s\pi } \over {2N}}} \over {1+4\alpha \sigma \sin ^{2}{{s\pi } \over {2N}}}}\;$

(8.29)

Rozwiązaniem całego zagadnienia (8.12) jest kombinacją liniową względem wskaźnika "s" (8.27):

$V_{i,j}=\sum _{s=1}^{N-1}a_{s}\lambda _{s}^{j}\sin {{s\pi i} \over {N}}\;$

(8.30)

Funkcja G_i, czyli funkcja (8.30) zapisana dla j=0, spełnia równanie zapisane bez udziału współczynników λ_s, ale przy udziale parametru N (która jest ilością punktów na osi iskowej w przedziale zmienności o długości "l" i "i", który jest numerem punktu na osi iksowej):

$G_{i}=V_{i,0}=\sum _{s=1}^{N-1}a_{s}\sin {{s\pi i} \over {N}}\;$

(8.31)

Mnożymy wyrażenie (8.29) przez $_{\sin {{r\pi i} \over {N}}}\;$ oraz z sumując tak otrzymane równanie:

$\sum _{i=1}^{N-1}G_{i}\sin {{r\pi i} \over {N}}=\sum _{i=1}^{N-1}\sum _{s=1}^{N-1}a_{s}\sin {{s\pi i} \over {N}}\sin {{r\pi i} \over {N}}\Rightarrow \sum _{i=1}^{N-1}G_{i}\sin {{r\pi i} \over {N}}=\sum _{s=1}^{N-1}a_{s}\sum _{i=1}^{N-1}\sin {{s\pi i} \over {N}}\sin {{r\pi i} \over {N}}\;$

(8.32)

Wyznaczmy teraz sumę, którą wykorzystamy w obliczeniach (8.32), do którego wykorzystamy obliczenia przeprowadzone w punkcie (2.49):

$\sum _{i=0}^{N-1}\sin {{s\pi i} \over {N}}\sin {{r\pi i} \over {N}}={{1} \over {2}}\sum _{i=0}^{N-1}\left[\cos {{(s-r)\pi i} \over {N}}-\cos {{(s+r)\pi i} \over {N}}\right]={{1} \over {2}}N\delta _{rs}\;$

(8.33)

Wykorzystując obliczenia przeprowadzone w punkcie (8.33) obliczenia (8.32) możemy dokończyć do końca i przepisać gotowy wynik na współczynnik a_s:

$a_{s}={{2} \over {N}}\sum _{i=1}^{N-1}G_{i}\sin {{r\pi i} \over {N}}{\mbox{ dla }}r=1,..,N-1\;$

(8.34)

Funkcja (8.28) jest określona przy pomocy λ_s (8.29) i przy pomocy współczynników a_s (8.34). Zbudujmy iloczyn skalarny oparty o iloczyn funkcji V i W przy kroku h:

$\langle V,W\rangle =\sum _{i=1}^{N-1}hV(x_{i})W(x_{i})\;$

(8.35)

Jest ona różnicowym odpowiednikiem $\langle V,W\rangle ={\sqrt {\int _{0}^{1}V(x)W(x)dx}}\;$. Weźmy sobie ciąg funkcji zależnych od s, który definiujemy w zależności od "l" przedstawiający długość osi iksowej, dla której jest słuszne rozwiązanie X^r(x):

$X^{s}(x)={\sqrt {{2} \over {l}}}\sin s\pi x{\mbox{ dla }}s=1,..,N-1\;$

(8.36)

Ciąg funkcji X¹(x), X²(x),...,X^N-1 jest ortogonalny na zbiorze punktów $x_{i}={{i} \over {N}}\;$ dla i=1,...,N-1, tzn.:

$\langle X^{n},X^{m}\rangle =\sum _{i=1}^{N-1}hX^{n}(x_{i})X^{m}(x_{i})={{2h} \over {l}}\sum _{i=1}^{N-1}\sin {{n\pi i} \over {N}}\sin {{m\pi i} \over {N}}=\;$
$={{2h} \over {l}}{{1} \over {2}}N\delta _{n,m}={{Nh} \over {l}}\delta _{n,m}=\delta _{n,m}\;$

(8.37)

Funkcję (8.28) zapisujemy przy pomocy funkcji (8.36) i względem wskaźników "j" i "i" oraz od długości przedziału zmienności iksów:

$V_{i,j}=\sum _{s=1}^{N-1}a_{s}\lambda _{s}^{j}\sin {{s\pi i} \over {N}}={\sqrt {{l} \over {2}}}\sum _{s=1}^{N-1}\lambda _{s}^{j}{\sqrt {{2} \over {l}}}\sin s\pi x_{i}={\sqrt {{l} \over {2}}}\sum _{s=1}^{N-1}\lambda _{s}^{j}X^{s}(x_{i})\;$

(8.38)

Napiszmy czemu jest równa norma ciągu (8.38), który napiszemy przy definicji iloczynu skalarnego (8.35):

$||V^{j}||_{L_{2}}^{2}=\langle V^{j},V^{j}\rangle ={{l} \over {2}}\sum _{s,k=1}^{N-1}a_{s}a_{k}\lambda _{s}^{j}\lambda _{k}^{j}\langle X^{i},X^{k}\rangle ={{l} \over {2}}\sum _{s=1}^{N-1}a_{s}^{2}\lambda _{s}^{2j}\;$

(8.39)

Gdy we wzorze (8.39) mamy j=0, to możemy napisać normę funkcji G wychodząc od normy funkcji V_i,j zależną od kwadratów współczynników a_s i od długości przedziału zmienności iksów, czyli "l":

$||G||_{L_{2}}^{2}={{l} \over {2}}\sum _{s=1}^{N-1}a_{s}^{2}\;$

(8.40)

Natomiast, gdy |λ_s|≤1, to normę ||V^j|| zapisujemy w zależności od współczynników a_s, jak się przekonamy, że opisywana norma jest zawsze mniejsza od normy ||G|| (8.40):

$||V^{j}||_{L_{2}}^{2}\leq {{l} \over {2}}\sum _{i=1}^{N-1}a_{s}^{2}=||G||_{L_{2}}^{2}\;$

(8.41)

Zbadajmy dla jakich współczynników σ,h,k zachodzi |λ_s|≤1, s=1,..,N-1, jeśli definicję λ_s zapisaliśmy w punkcie (8.27):

$-1\leq \lambda \leq 1\Rightarrow -1\leq 1-{{4\alpha \sin ^{2}{{s\pi } \over {2N}}} \over {1+4\alpha \sigma \sin ^{2}{{s\pi } \over {2N}}}}\leq 1\Rightarrow -1\leq 1-{{1} \over {{{1} \over {4\alpha \sin ^{2}{{s\pi } \over {2N}}}}+\sigma }}\leq 1\;$
$-{{1} \over {4\alpha \sin ^{2}{{s\pi } \over {2N}}}}-\sigma \leq {{1} \over {4\alpha \sin ^{2}{{s\pi } \over {2N}}}}+\sigma -1\leq {{1} \over {4\alpha \sin ^{2}{{s\pi } \over {2N}}}}+\sigma \Rightarrow \;$
$\Rightarrow 1-{{1} \over {2\alpha \sin ^{2}{{s\pi } \over {2N}}}}\leq 2\sigma \Rightarrow \sigma \geq {{1} \over {2}}-{{1} \over {4\alpha \sin ^{2}{{s\pi } \over {2N}}}}\;$

(8.42)

Z nierówności (8.40) możemy wywnioskować dla s=1,...,N-1 i z definicji funkcji sinus:

$\sigma \geq {{1} \over {2}}-{{1} \over {4\alpha }}={{1} \over {2}}-{{h^{2}} \over {4k}}\;$

(8.43)

Jeśli natomiast σ<1/2, to na podstawie obliczeń (8.43) ograniczenie na stosunek liczby k przez kwadrat liczby h²:

$4\sigma \geq 2-{{1} \over {\alpha }}\Rightarrow {{1} \over {\alpha }}\geq 2(1-2\sigma )\Rightarrow \alpha \leq {{1} \over {2(1-\sigma )}}\Rightarrow {{k} \over {h^{2}}}\leq {{1} \over {2(1-2\sigma )}}\;$

(8.44)

dla tego σ opisywana metoda jest stabilna względem warunku początkowego. Dla σ=0 okreslenie równania różnicowego jest zadaniem bardzo prostym. Już dla σ≥1/2 metoda różnicowa jest słuszna dla dowolnych h i k. Wykażemy, jeśli $_{\sigma \leq {{1} \over {2}}-{{1} \over {4\alpha }}-\epsilon }\;$ (gdzie ε jest dowolną stałą dodatnią), to ona jest niestabilna względem normy ||⋅||_L2. Obierzmy N na tyle duże, by spełniony warunek s₀<N, takie by był spełniony warunek |λ_s0|>1+ε₁, wtedy na podstawie (8.39) powiemy:

$||(V^{s_{0}})^{j}||_{L_{2}}^{2}=\lambda _{s_{0}}^{2j}{{l} \over {2}}\geq {{l} \over {2}}(1+\epsilon _{1})^{2j}\rightarrow \infty \;$

(8.45)

Gdy k=0, wtedy "j" jest bardzo duże, i dlatego wyrażenie (8.45) dąży do nieskończoności. Rozpatrzmy przypadek σ=0, wtedy drugą pochodną w równaniu przewodnictwa cieplnego $_{{\partial ^{2}u} \over {\partial u^{2}}}\;$ możemy przybliżyć wyrażeniem w postaci ilorazu różnicowego przy oznaczeniu U_i(t)=u(x,t):

${{1} \over {h^{2}}}\left[U_{i+1}(t)-2U_{i}(t)+U_{i-1}(t)\right]\;$

(8.46)

Równanie (8.1) w takim bądź razem zapisujemy przy rozbiciu drugiej pochodnej wielkości U(x,t) względem położenia na wyrazenie różnicowe (8.44):

$U_{i}^{'}(t)={{1} \over {h^{2}}}\left[U_{i+1}(t)-2U_{i}(t)+U_{i-1}(t)\right]\;$

(8.47)

Dla którego mamy warunki brzegowe U₀(t)=φ(t) i U_N(t)=ψ(t) dla 0≤t≤T i z warunkiem poczatkowym U_i(0)=G(x_I) dla i=1,...,N-1. Równanie (8.45) możemy zapisać w postaci macierzowej wzgledem wektora pionowego opisującego wielkość, którego elementami są U_i:

$\mathbf {U} ^{'}(t)={{1} \over {h^{2}}}\mathbf {M} \mathbf {U} (t)+\mathbf {F} (t)\;$

(8.48)

z warunkiem początkowym U(0)=G. Macierze występujące w (8.48) są to obiekty zdefiniowane jako:

$\mathbf {M} ={\begin{bmatrix}-2&1&&&0&\\1&-2&1&&\\&\ddots &\ddots &\ddots \\&0&1&-2&1\\&&&1&-2\\\end{bmatrix}}\;$

(8.49)

$\mathbf {U} (t)={\begin{bmatrix}U_{1}(t)\\U_{2}(t)\\\vdots \\U_{N-1}\end{bmatrix}}\;$

(8.50)

$\mathbf {F} (t)={\begin{bmatrix}f(x_{1},t)+{{1} \over {h^{2}}}\phi (t)\\f(x_{2},t)\\\vdots \\f(x_{N-2},t)\\f(x_{N-1},t)+{{1} \over {h^{2}}}\psi (t)\\\end{bmatrix}}\;$

(8.51)

$\mathbf {G} ={\begin{bmatrix}g(x_{1})\\g(x_{2})\\\vdots \\g(x_{N-1})\end{bmatrix}}\;$

(8.52)

Równania różnicowe i ich zbieżność względem parametrów siatki względem osi położenia i czasu

Weźmy sobie definicję $F_{i,j+1/2}=L_{k,h}^{\sigma }(v_{ij}-u_{ij})\;$, gdy jest spełnione G_i=0, wtedy możemy odjąć równanie (8.12) i (8.13) od równania (8.6), w ten sposób możemy otrzymać ograniczenie od góry na |F_i,j+1/2| dla odpowiednich σ przy użyciu stałych c₁, c₂:

$|F_{i,j+1/2}|\leq \left\{{\begin{matrix}c_{1}(h^{2}+k^{2})&\sigma =1/2&u\in C_{3}^{4}({\overline {\mathbf {D} }})\\c_{2}(h^{2}+k)&\sigma \neq 1/2&u\in C_{2}^{4}({\overline {\mathbf {D} }})\end{matrix}}\right.\;$

(8.53)

Patrząc na udowodnione twierdzenie (8.26) otrzymane w końcowym wniosku w tym punkcie, a także wykorzystując (8.53), otrzymujemy bardzo ważny wniosek na oszacowania różnicy v_ij-u_ij względem normy ||⋅||_c:

$\max _{(i,j)\in {\hat {\mathbf {\omega } }}}|v_{ij}-u_{ij}|=\max _{j}||v^{j}-u^{j}||_{c}\leq T\max _{0\leq s\leq j-1}||F^{s+1/2}||\leq T\max _{(i,j)\in {\hat {\mathbf {\omega } }}}F_{i,j+1/2}\leq \;$
$\leq \left\{{\begin{matrix}c_{1}T(h^{2}+k^{2})&\sigma =1/2&u\in C_{3}^{4}({\overline {\mathbf {D} }})\\c_{2}T(h^{2}+k)&\sigma \neq 1/2&u\in C_{2}^{4}({\overline {\mathbf {D} }})\end{matrix}}\right.\;$

(8.54)

Norma błędu $L_{h,k}^{\sigma }(v_{ij}-u_{ij})\;$dąży do zera, gdy wielkości "h" i "k" dążą do zera, przy zastąpieniu równania różniczkowego równaniem różnicowym, co na podstawie (8.54) możemy powiedzieć:

$\max _{j}||v^{j}-u^{j}||{\xrightarrow[{h,k\rightarrow 0}]{}}0\;$

(8.55)

Zbieżność przy nieskończonym zagęszczaniu siatki

Równanie różnicowe (8.7) podobnie się rozwiązuje jak równanie (8.14), to na podstawie (8.30) i wzoru na λ_s (8.29) możemy otrzymać wzór na v_ij przy założeniu $\alpha ={{h} \over {h^{2}}}=\operatorname {const} \;$ i σ=0, wiedząc przy tym, że x=ih i t=jk, dalej przechodząc do granicy h→0 i k→0, otrzymujemy ważne dokładne rozwiązanie:

$v_{ij}=\left(1-4\alpha \sin ^{2}{{\pi } \over {2N}}\right)^{i}\sin {{\pi i} \over {N}}\Rightarrow \lim _{\overset {h\rightarrow 0}{\underset {k\rightarrow 0}{}}}v_{ij}=\lim _{\overset {h\rightarrow 0}{\underset {k\rightarrow 0}{}}}\left(1-4\alpha \sin ^{2}{{\pi {\sqrt {k}}} \over {2{\sqrt {\alpha }}}}\right)^{{t} \over {k}}\sin \pi x=\;$
$=\lim _{\overset {h\rightarrow 0}{\underset {k\rightarrow 0}{}}}(1-\pi ^{2}k)^{{t\pi ^{2}} \over {\pi ^{2}k}}\sin \pi x=e^{-\pi ^{2}t}\sin \pi x\;$

(8.56)

Rozwiązanie równania (8.1) przy f(x,t)=0 jest u(x,t)=e^-π2tsinπx.

Rozwiązywanie różnicowego równania według schematu blokowego

Równanie (8.7) dzielimy przez -ασ i w rezultacie otrzymujemy następne równanie przy definicji $s={{h^{2}} \over {\sigma k}}\;$:

$v_{i+1,j}-(2+s)v_{i,j}+v_{i-1,j}=-r_{i,j}\;$
$r_{i,j}={{1-\sigma } \over {\sigma }}v_{i+1,j-1}+\left[s-2{{1-\sigma } \over {\sigma }}\right]v_{i,j-1}+{{1-\sigma } \over {\sigma }}v_{i-1,j-1}-{{h^{2}} \over {\sigma }}f_{i,j-1/2}\;$
$v_{0j}=\varphi _{j}{\mbox{, }}v_{N,j}=\psi _{j}\;$

(8.57)

Rozwiązanie (8.57) będziemy wyznaczać kolejno warstwami dla t=jk. Dla każdej warstwy mamy do rozwiązania N-1 równań tworzące jeden układ równań. Tutaj będziemy rozważać "j" jako ustalona liczba. Ponieważ równanie (8.57) zapisujemy względem macierzy trójdiagonalnej, który zapisujemy w formie (5.93). Podczas działania macierzy trójdiagonalnej na macierz o elementach v_i,j (które będziemy wyznaczać) otrzymujemy wektor pionowy o elementach r_ij. Z tożsamości (5.103) otrzymujemy zależność łączącą v_i,j z v_i+1,j poprzez parametry p_i i t_i,j:

$v_{i,j}=p_{i}v_{i+1,j}+t_{i,j}{\mbox{ dla }}i=N-1,..,1\;$

(8.58)

wyznaczmy elementy p_i wykorzystując wzór (5.98) oraz (5.99), wtedy możemy napisać związek rekurencyjny p_i względem p_i-1 w zależności dodatkowo od parametru "s" zdefiniowanego powyżej:

$p_{i}=-{{c_{i}} \over {u_{i}}}=-{{1} \over {u_{i}}}=-{{1} \over {-(2+s)-l_{i}c_{i-1}}}=-{{1} \over {-(2+s)-{{1} \over {u_{i-1}}}}}={{1} \over {2+s-p_{i-1}}}\;$

(8.59)

$p_{0}=0\;$

(8.60)

Dalej jeśli wykorzystamy równanie (5.101), to otrzymamy zależność rekurencyjną t_i,j w zależności od t_i-1,j poprzez parametr r_i,j zdefiniowanej w punkcie (8.57) w jego drugiej linijce:

$t_{i,j}={{y_{i}} \over {u_{i}}}=-p_{i}y_{i}=-p_{i}\left(-r_{ij}-{{y_{i-1}} \over {u_{i-1}}}\right)=p_{i}(t_{i-1,j}+r_{i,j})\;$

(8.61)

Patrząc na (8.57) na trzecią linijkę (pierwszy wzór w tej linijce), a także na wzór (8.60) i (8.58), dochodzimy natychmiast, że:

$t_{0,j}=\varphi _{j}\;$

(8.62)

Drgająca struna

Równanie struny jednorodnej jako przykład równania hiperbolicznego przedstawimy jako równanie różniczkowe, w której wyznaczaną funkcją jest u(x,t), którą wyznaczamy względem funkcji f(x,t):

$Lu={{\partial ^{2}u} \over {\partial t^{2}}}-{{\partial ^{2}u} \over {\partial x^{2}}}=f(x,t)\;$

(8.63)

Równanie (8.63) jest określone w dziedzinie napisanej wewnątrz prostokąta ${\overline {\mathbf {D} }}=\{(x,t):0\leq x\leq l;,0\leq t\leq T\}\;$. Warunkiem początkowym dla rozwiązania struny jednorodnej (8.63) jest u(x,0)=g(x) i ${{\partial u(x,0)} \over {\partial t}}=F(x)\;$ dla 0≤x≤l, a warunkami brzegowymi są u(0,t)=φ(t), u(l,t)=ψ(t) dla 0≤t≤T. Zdefiniujmy operator $L_{h,k}^{\sigma }\;$, który działa na funkcję v_ij wedle:

$L_{h,k}^{\sigma }v_{ij}={{v_{i,j+1}-2v_{i,j}+v_{i,j-1}} \over {k^{2}}}-\mathbf {\Lambda } \left(\sigma v_{i,j+1}+(1-2\sigma )v_{i,j}+\sigma v_{i,j-1}\right)=f_{i,j}\;$

(8.64)

Zakładamy, że naszymi warunkami brzegowymi są v_i,0=g_i dla i=1,...,N-1 i v_0,j=φ_j, v_N,j=ψ_j dla j=0,1,...,M, gdzie g_i=g(x_i), φ_i=φ(t_j) oraz ψ_j=ψ(t_j). Jeśli obierzemy funkcję $u\in C_{4}^{4}({\overline {\mathbf {D} }})\;$ dla dowolnego parametru σ, która nie zależy od parametru podziału siatki, tzn. od h, k, wtedy możemy napisać:

$L_{h,k}^{\sigma }u_{i,j}=f_{i,j}+O(h^{2}+k^{2})\;$

(8.65)

Rozważmy przypadek brzegowy dla czasu równego zero, co na tej podstawie równania (8.63) iloraz różnicowy funkcji u(x,t) piszemy:

${{u_{i,1}-u_{i,0}} \over {k}}={{\partial u_{i,0}} \over {\partial t}}+{{k} \over {2}}{{\partial ^{2}u_{i,0}} \over {\partial t^{2}}}+O(k^{2})={{\partial u_{i,0}} \over {\partial t}}+{{k} \over {2}}\left(f_{i,0}+{{\partial u_{i,0}} \over {\partial x^{2}}}\right)=\;$
$={{\partial u_{i,0}} \over {\partial t}}+{{k} \over {2}}\left(f_{i,0}+g_{i}^{''}\right)+O(k^{2})=F_{i}+{{k} \over {2}}\left(f_{i,0}+g_{i}^{''}\right)+O(k^{2})\;$

(8.66)

Jeśli będziemy rozpatrywać funkcje v_i,j dla warunku początkowego, to na podstawie (8.66) możemy napisać równość:

${{v_{i,1}-v_{i,0}} \over {k}}={{\partial v_{i,0}} \over {\partial t}}+{{k} \over {2}}\left(f_{i,0}+g_{i}^{''}\right)=F_{i}+{{k} \over {2}}\left(f_{i,0}+g_{i}^{''}\right)\;$

(8.67)

Równanie (8.64) będziemy rozwiązywać warstwami dla i=1,..,N-1 i go możemy przepisać w postaci dwóch równości przy pewnej definicji r_ij:

$\left({{1} \over {k^{2}}}-\sigma \mathbf {\Lambda } \right)v_{i,j+1}=r_{ij}\;$

(8.68)

$r_{ij}=\left({{2} \over {k^{2}}}+(1-2\sigma )\mathbf {\Lambda } \right)v_{i,j}+\left(-{{1} \over {k^{2}}}+\sigma \mathbf {\Lambda } \right)v_{i,j-1}+f_{ij}\;$

(8.69)

Rozwiązywanie równań Poissona metodami różnicowymi

Weźmy najprostszy przykład równania Poissona obowiązującego w przestrzeni trójwymiarowej:

$Lu={{\partial ^{2}u} \over {\partial x^{2}}}+{{\partial ^{2}u} \over {\partial y^{2}}}=f(x,y)\;$

(8.70)

Funkcja f(x,y) w powyższym równaniu powinna być ciągła na obszarze ${\overline {\mathbf {D} }}\;$. Rozpatrzmy dla równania (8.69) zagadnienie brzegowe pierwszego rodzaju (zagadnienie Dirichleta), który wewnątrz tego obszaru u(x,y) spełnia równanie Poissona, a na brzegu tego obszaru jest spełnione zagadnienie $u|_{\mathbf {\Gamma } }=\varphi \;$, gdzie funkcja φ jest określona na brzego Γ. Rozpatrzmy drugie pochodne na obszarze Γ.

(Rys. 8.1) Siatka prostokątna a równanie Poissona

Rozpatrzmy drugie pochodne w węźle zerowym, tzn. ${{\partial ^{2}v} \over {\partial t^{2}}}\;$, ${{\partial ^{2}v} \over {\partial y^{2}}}\;$, które zapisujemy sposobem:

${{\partial ^{2}v_{0}} \over {\partial t^{2}}}={{v_{1}-2v_{0}+v_{3}} \over {h^{2}}}\;$

(8.71)

${{\partial ^{2}v_{0}} \over {\partial y^{2}}}={{v_{2}-2v_{0}+v_{4}} \over {k^{2}}}\;$

(8.72)

Całe równanie Poissona (8.69) dla węzła zerowego możemy zapisać jako:

$L_{h,k}v_{0}={{v_{1}-2v_{0}+v_{3}} \over {h^{2}}}+{{v_{2}-2v_{0}+v_{4}} \over {k^{2}}}=f_{0}\;$

(8.73)

Jesli natomiast będziemy rozważać funkcje $u\in C^{4}({\overline {\mathbf {D} }})\;$, gdzie: $C^{m}({\overline {\mathbf {D} }})\;$, jest klasą funkcji mające ciągłe pochodne do rzęd m-tego włącznie w zbiorze ${\overline {\mathbf {D} }}\;$, to możemy napisać:

$L_{h,k}u_{0}=f_{0}+O(h^{2}+k^{2})\;$

(8.74)

Użyjmy dziewięć punktów opisane na rysunku obok, dzięki któremu równanie aproksymacyjne, które po odpowiednim przegrupowaniu wyrazów, możemy napisać:

${{5} \over {3}}\left({{1} \over {h^{2}}}+{{1} \over {k^{2}}}\right)v_{0}={{1} \over {6}}\left({{5} \over {h^{2}}}-{{1} \over {k^{2}}}\right)(v_{1}+v_{3})+{{1} \over {6}}\left({{5} \over {k^{2}}}+{{1} \over {h^{2}}}\right)(v_{2}+v_{4})+\;$
$+{{1} \over {12}}\left({{1} \over {h^{2}}}+{{1} \over {k^{2}}}\right)(v_{5}+v_{6}+v_{7}+v_{8})+f_{0}+{{h^{2}} \over {12}}(f_{1}-2f_{0}+f_{3})+{{k^{2}} \over {12}}\left(f_{2}-2f_{0}+f_{4}\right)\;$

(8.75)

(Rys. 8.2) Siatka trójkątna

(Rys. 8.3) Siatka trójkątna-wycinek

Rozpatrzmy siatkę trójkątną wyglądającą jak w powyższych rysunkach, wtedy równanie różnicowe symbolizujące równanie Poissona (8.69) jest:

${{2} \over {3h^{2}}}(v_{1}+v_{2}+v_{3}+v_{4}+v_{5}+v_{6}-6v_{0})=f_{0}\;$

(8.76)

Problem Dirichleta w aproksymacji różnicowej

(Rys. 8.4) Krzywa z zaznaczeniem punktów należących do krzywej Γ i punktów wewnętrznych D

Określmy sobie zbiór $\mathbf {D} \cup \mathbf {\Gamma } \;$ składających się z punktów wewnętrznych D i punktów krzywej Γ, wtedy możemy poszukiwać rozwiązania zagadnienia Dirichleta, pokrywając płaszczyznę prostokątną siatką prostych:

$x=a+ih{\mbox{, }}y=b+jh\;$

(8.77)

Określmy przez ${\overline {\mathbf {\omega } }}\;$ zbiór węzłów należących do zbioru wyżej określonego, które dzielimy na zbiór węzłów wewnętrznych ω i na zbiór węzłów brzegowych. Dwa punkty nazywamy sąsiednimi, gdy jest spełniony warunek:

$\sum _{s=1}^{2}|i_{s}-j_{s}|=1\;$

(8.78)

Węzeł nazywamy wewnętrzny, jeśli każdy jego sąsiad leży w zbiorze ${\overline {\mathbf {\omega } }}\;$, a pozostałe węzły nazywamy brzegowymi. Zbiór węzłów brzegowych podzielmy na dwa zbiory, zbiór γ₁, które leżą bezpośrednio na krzywej Γ, a pozostałe węzły brzegowe oznaczmy przez γ₂. Przez zamalowane kółeczko na rysunku obok będziemy oznaczać węzły ω, a przez niezamalowane kółeczko węzły γ₁, a przez niezmalowaną gwiazdkę węzły γ₂. Jeśli punkt $A\in \mathbf {\gamma } _{1}\;$, to na pewno powiemy:

$v(A)=\phi (A)\;$

(8.79)

(Rys. 8.5) Interpolacja wartości

Dla punktów węzłowych siatki w obszarze D równanie różnicowe symbolizujące równanie Poissona piszemy:

$L_{h,k}v_{ij}={{v_{i+1,j}-2v_{ij}+v_{i-1,j}} \over {h^{2}}}+\;$
$+{{v_{i,j+1}-2v_{ij}+v_{i,j-1}} \over {k^{2}}}=f_{ij}\;$

(8.80)

Dla punktów brzegowych, ale nieleżących na krzywej Γ, czyli należącej do obszaru γ₂ można w przybliżeniu napisać v_ij=φ(B_ij). Dla punktów brzegowych można policzyć wartość v_ij znając wartość funkcji v(x,y) na krzywej Γ (punkt B) i znając wartość dla punktów wewnętrznych (punkt C) według interpolacji liniowej jak to podał Colatz z błędem O(h²+k²) według wzoru:

$v(A)={{1} \over {h+h_{1}}}\left(h_{1}v(C)+h\phi (B)\right)\;$

(8.81)

(Rys. 8.6) Zastosowanie równania różnicowego dla poszczególnych punktów

Weźmy sobie punkt A, który jest punktem brzegowym należącym do zbioru γ₂. Weźmy dodatkowo cztery punkty, które otaczają ten punkt A, które z tych punktów są albo węzłami, albo punktami przecięcia siatki z krzywą Γ, wtedy równanie różnicowe możemy zapisać przez:

$l_{h,k}v(A)={{1} \over {h}}\left({{v(B)-v(A)} \over {h_{1}}}-{{v(A)-v(D)} \over {h_{2}}}\right)+\;$
$+{{1} \over {k}}\left({{v(C)-v(A)} \over {k_{1}}}-{{v(A)-v(E)} \over {k_{2}}}\right)=f(A)\;$

(8.82)

We wzorze (8.82) w punktach, które pokrywają się z punktami na krzywej γ podstawiamy wartości funkcji φ w tychże punktach. Równanie (8.82) możemy zapisać w sposób bardziej przekształcony pisząc je:

$v(A)\left[{{1} \over {h_{1}h}}+{{1} \over {h_{2}h}}+{{1} \over {kk_{1}}}+{{1} \over {kk_{2}}}\right]-{{1} \over {h_{1}h}}v(B)-{{1} \over {hh_{2}}}v(D)-{{1} \over {kk_{1}}}v(c)-{{1} \over {kk_{2}}}v(E)=f(A)\;$

(8.83)

Rozwiązywanie równań różnicowych w wykorzystaniem zbieżności jednostajnej

Jeśli będziemy korzystać z równości (8.79) i v_ij=φ(B_ij), wtedy oszacowanie maksymalne pomiędzy wartością u_ij, a v_ij jest:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|u_{ij}-v_{ij}|\leq N_{1}(h+k)\;$

(8.84)

a jeśli weżniemy pod uwagę (8.79), (8.81) i (8.79), (8.82), wtedy oszacowanie błędu pomiędzy wartościami u_ij i v_ij jest:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|u_{ij}-v_{ij}|\leq N_{2}(h^{2}+k^{2})\;$

(8.85)

Dla węzłów, których wartości są v_ij możemy zbudować twierdzenia, które jest bardzo łatwo udowodnić na podstawie rozdziału Stabilność rozwiązań równań cząstkowych z warunkami brzegowymi i jego zbieżność, które tutaj przedstawimy bez dowodu:

Twierdzenie pierwsze
Dla v niestałego wyrażenie L_h,k≥0 (L_h,k≤0) funkcja w węzłach wewnętrznych nie może przyjmować wartości największej (najmniejszej).

Twierdzenie drugie
Jeśli są spełnione warunki na węzłów wewnętrznych |L_h,kv_ij|≤-L_h,kY_ij i |v_ij|≤Y_ij dla węzłów brzegowych, to wtedy na pewno dla węzłów brzegowych i wewnętrznych spełniony jest warunek |v_ij|≤Y_ij.

Równanie (8.79) ma jednoznaczne rozwiązanie w funkcjach v_ij, w tym celu należy rozwiązać równanie jednorodne:

$L_{h,k}v_{ij}=0\;$

(8.86)

i stwierdzić czy ono ma tylko rozwiązania v_ij zerowe, ale przy wykorzystaniu twierdzenia pierwszego stwierdzamy, że nie ma ani rozwiązań dodatnich, ani ujemnych, zatem jedynymi rozwiązaniami są rozwiązania zerowe. Objeżmy sobie funkcję v=z+q, wtedy możemy powiedzieć dla funkcji "z", że L_h,kz_i,j=0 dla punktów wewnętrznych i z_ij=φ(B_ij) dla punktów brzegowych, a także dla funkcji "q" zachodzą dwa warunki L_h,kq_ij=f_ij dla punktów wewnętrznych i q_ij=0 dla punktów brzegowych. Na podstawie powyższych wniosków możemy powiedzieć:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|z_{ij}|\leq \max _{(i,j)\in \mathbf {\omega } }|z_{ij}|\leq \max _{A\in \mathbf {\Gamma } }\phi (A)\;$

(8.87)

Przyjmijmy teraz funkcję Y(x,y)=K(r²-x²-y²), gdzie przez R będziemy oznaczać promieniem koła, w której zawarty jest obszar D, wtedy na podstawie tego napiszemy L_h,kY_ij=-4K. Przyjmijmy oznaczenie:

$K=\max _{(i,j)\in \mathbf {\omega } }{{|f_{ij}|} \over {4}}\;$

(8.88)

wtedy możemy przyjąć oznaczenie na q_ij na podstawie twierdzenia drugiego:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|q_{ij}|\leq \max _{(i,j)\in {\overline {\mathbf {\omega } }}}|Y_{ij}|\leq R^{2}\max _{(i,j)\in {\overline {\mathbf {\omega } }}}{{f_{ij}} \over {4}}={{R^{2}} \over {4}}\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|f_{ij}|\leq {{R^{2}} \over {4}}\max _{A\in D}|f(A)|\;$

(8.89)

zatem funkcję v_ij według definicji samego v podanego powyżej możemy napisać:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|v_{ij}|\leq \max _{A\in \mathbf {\Gamma } }|\phi (A)|+{{R^{2}} \over {4}}\max _{A\in D}|f(A)|\;$

(8.90)

Patrząc na (8.89) udowodniliśmy, że wyrazy v_ij zależą w sposób ciągły od funkcji f i φ. Weźmy sobie funkcję błędu ε=u-v, wtedy napiszemy L_h,kε_ij=F_ij dla punktów wewnętrznych i ε_ij=ψ_ij dla punktów brzegowych, wtedy oszacowania na funkcje F_ij i ψ_ij możemy napisać przy definicji stałej M_s równą $_{M_{s}=\max \left\{\max _{(x,y)\in \mathbf {D} }\left|{{\partial ^{s}u} \over {\partial x^{s}}}\right|,\max _{(x,y)\in \mathbf {D} }\left|{{\partial ^{s}u} \over {\partial y^{s}}}\right|\right\}}\;$:

$|F_{ij}|\leq {{M_{4}} \over {12}}(h^{2}+k^{2})\;$ (8.91)

$|\psi _{ij}|\leq M_{1}(k+h)\;$ (8.92)

wtedy patrząc na wzór (8.90), a także na podstawie (8.91) i (8.92) możemy napisać na maksymalną wartość bezwzględną różnicy pomiędzy u_ij i v_ij według:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|u_{ij}-v_{ij}|\leq M_{1}(k+h)+{{R^{2}} \over {48}}M_{4}(h^{2}+k^{2})\;$

(8.93)

Jeśli brzeg jest prostokątem, to można tak wybrać brzeg tego prostokąta by on leżał na liniach tworzących naszą prostokątną siatkę, tzn. ψ_ij=0, to wtedy dla tego przypadku (8.93) możemy przepisać dla przypadku M₁ równą zero. Niech kolejno $v,{\tilde {v}}\;$ będą rozwiązaniami równania różnicowego przy $f,\phi ,{\tilde {f}},{\tilde {\phi }}\;$, wtedy zaburzenie wielkości v_ij na podstawie (8.90) piszemy wzorem:

$\max _{(i,j)\in {\overline {\mathbf {\omega } }}}|v_{ij}-{\tilde {v}}_{ij}|\leq \max _{A\in \mathbf {\Gamma } }|\phi (A)-{\tilde {\phi }}(A)|+{{R^{2}} \over {4}}\max _{A\in D}|f(A)-f({\tilde {A}})|\;$

(8.94)

Rozwiązywanie równania przewodnictwa cieplnego metodą prostych

Rozpatrzmy równanie przewodnictwa cieplnego, którego rozwiązaniem są funkcje wektorowe u=u(x,t) spełnionego wewnątrz prostokąta ${\overline {\mathbf {D} }}=\{(x,t):0\leq x\leq l,0\leq t\leq T\}\;$, zapisaną w postaci:

$L\mathbf {u} ={{\partial \mathbf {u} } \over {\partial t}}-{{\partial ^{2}\mathbf {u} } \over {\partial x^{2}}}=f(x,t)\;$

(8.95)

z warunkiem początkowym u(x,0)=g₁(x) dla 0≤x≤t i z warunkami brzegowymi u(0,t)=φ(t) i u(l,t)=ψ(t) dla 0≤t≤T. Rozwiązania zagadnienia będziemy poszukiwali metodą prostych x=x_i=ih przy zdefiniowanym h=l/N, gdzie i=0,1,..,N. Drugą pochodną wielkości u względem położenia w równaniu (8.95) zastępujemy wyrażeniem różnicowym przy pomocy prostych x=x_i dla i=1,...,N-1, wtedy po tej operacji piszemy:

$L_{h,k}v_{i}(t)=v^{'}(t)-{{1} \over {h^{2}}}\left[v_{i+1}(t)-2v_{i}(t)+v_{i-1}(t)\right]=f_{i}(t)\;$

(8.96)

przy warunkach początkowych v_i(0)=g₁(x_i) i przy warunkach brzegowych v₀(t)=φ(t), v_N(t)=ψ(t) dla 0≤t≤T. Równanie (8.96) możemy zapisać w postaci wektorowej wraz z warunkiem początkowym:

$\mathbf {V} ^{'}(t)={{1} \over {h^{2}}}\mathbf {A} \mathbf {V} (t)+\mathbf {F} (t){\mbox{ gdzie }}\mathbf {V} (0)=\mathbf {G} _{1}\;$

(8.97)

Metoda prostych polega na rozwiązaniu zagadnienia rozpisanego powyższym równaniem w postaci składowych V(t). Rozwiązanie równania macierzowego (8.97) jest zbieżne do rozwiązania dokładnego z dokładnością O(h²) i macierze tam występujące piszemy wzorami:

$\mathbf {A} ={\begin{bmatrix}-2&1&&0\\1&-2&1\\&\ddots &\ddots &\ddots \\&&-2&1\\0&&1&-2\end{bmatrix}}\;$

(8.98)

$\mathbf {V} (t)={\begin{bmatrix}v_{1}(t)\\v_{2}(t)\\\vdots \\v_{N-2}(t)\\v_{N-1}(t)\end{bmatrix}}\;$

(8.99)

$\mathbf {F} (t)={\begin{bmatrix}f_{1}(t)+{{1} \over {h^{2}}}\phi (t)\\f_{2}(t)\\\vdots \\f_{N-2}(t)\\f_{N-1}(t)+{{1} \over {h^{2}}}\psi (t)\end{bmatrix}}\;$

(8.100)

$\mathbf {G} _{1}={\begin{bmatrix}g_{1}(x_{1})\\g_{1}(x_{2})\\\vdots \\g_{1}(x_{N-2})\\g_{1}(x_{N-1})\end{bmatrix}}\;$

(8.101)

Jeśli mamy rozwiązanie równania różniczkowego (8.95) w postaci składowych o elementach u(x_i,t) i rozwiązania równania różnicowego (8.97) o elementach v_i(t), to wtedy dla $\mathbf {u} \in C_{2}^{4}({\overline {\mathbf {D} }})\;$:

$|u(x_{i},t)-v_{i}(t)|\leq Mh^{2}\;$

(8.102)

Rozwiązywanie równania drgania struny metodą prostych

Rozpatrzmy równanie drgań struny, którego rozwiązaniem są funkcje wektorowe u=u(x,t) spełnionego wewnątrz prostokąta $_{{\overline {\mathbf {D} }}=\{(x,t):0\leq x\leq l,0\leq t\leq T\}}\;$, zapisaną w postaci:

$L\mathbf {u} ={{\partial ^{2}\mathbf {u} } \over {\partial t^{2}}}-{{\partial ^{2}\mathbf {u} } \over {\partial x^{2}}}=f(x,y)\;$

(8.103)

z warunkami początkowymi u(x,0)=g₁(x), ${{\partial \mathbf {u} (x,0)} \over {\partial t}}=g_{2}(x)\;$ dla 0≤x≤l i z warunkami brzegowymi u(0,t)=φ(t), u(l,t)=φ(t) dla 0≤t≤T. Równanie różniczkowe (8.103) zastępujemy równaniem różnicowym o skoku iksowym h:

$L_{h}v_{i}(t)=v_{i}^{''}(t)-{{1} \over {h^{2}}}\left[v_{i+1}(t)-2v_{i}(t)+v_{i-1}(t)\right]=f_{i}(t)\;$

(8.104)

przy warunkach brzegowych v_i(0)=g₁(x_i), v'_i(0)=G₂(x_i). Równanie różnicowe (8.104) możemy zastąpić równaniem wektorowym zapisanym w postaci jednego równania:

$V^{''}(t)={{1} \over {h^{2}}}\mathbf {A} V(t)+\mathbf {F} (t){\mbox{ przy }}\mathbf {V} (0)=\mathbf {G} _{1}{\mbox{ i }}\mathbf {V} ^{'}(0)=\mathbf {G} _{2}\;$

(8.105)

Wektory występujące w równaniu (8.105) piszemy takimi samymi wzorami co w punkcie (8.98), (8.99), (8.100) , a wektor G₂ podobnie definiujemy jak wektor (8.99), ale ze wskaźnikiem dwa zamiast jeden. Jeśli mamy rozwiązanie równania różniczkowego (8.103) w postaci składowych o elementach u(x_i,t) i rozwiązania równania różnicowego (8.104) o elementach v_i(t), to wtedy dla $u\in C_{4}^{4}({\overline {\mathbf {D} }})\;$:

$|u(x_{i},t)-v_{i}(t)|\leq M_{1}h^{2}\;$

(8.106)