yj+1 = yj + hf (tj, yj)
y0 = x0
Wzór ten opisuje rekurencyjnie ciąg punktów yi, który ma być ciągiem przybliżeń wartości rzeczywistego rozwiązania w punktach ti. Podamy warunek, jaki musi spełniać zagadnienie, aby przy zwiększeniu liczby tych punktów, czyli zmniejsza-niu długości przedziału h, przybliżenie było coraz lepsze.
4. Z otrzymanego ciągu punktów budujemy funkcję y, która będzie przybliżeniem całego rozwiązania. Można opisać ją wzorem:
(
y(tj) = yj
y(t) = yj + (t − tj)yj+1−yj
h = yj+ (t − tj)f (tj, yj) dla t ∈ [tj, tj+1] Wykresem takiej funkcji jest łamana, zwana łamaną Eulera.
W dalszych rozważaniach zakładamy, że metoda Eulera daje nam „od razu” przy-bliżenie rozwiązania na całym przedziale [a, b], a nie na jego części, tak jak to było po-dane na początku tego rozdziału. Wiemy że da się otrzymać takie przybliżenie na całym przedziale, więc uproszczenie to nie wpływa na poprawność formułowanych twierdzeń i faktów.
8.3 Zbieżność metody Eulera
Definicja 8.3.1. Mówimy, że metoda Eulera jest zbieżna, jeśli:
lim
N →∞ max
j=0,...,N|yj− x(tj)| = 0.
Gdzie yj jest punktem zdefiniowanym zgodnie z metodą Eulera, natomiast x(tj) jest wartością rzeczywistego rozwiązania w punkcie tj.
Sens tej definicji jest taki: metoda Eulera jest zbieżna, jeśli przy zwiększeniu liczby punktów tj różnica między prawdziwym rozwiązaniem, a rozwiązaniem przybliżonym maleje. Czyli im podział będziesz gęstszy, tym przybliżenie lepsze.
Twierdzenie 8.3.2 (o zbieżności metody Eulera). Załóżmy, że funkcja f ∈ C([a, b] ×
R, R), oraz istnieje takie L > 0, że ∀t∈[a,b]∀x,y∈R|f (t, x) − f (t, y)| ¬ L|x − y|. Wówczas metoda Eulera dla zagadnienia (P C) jest zbieżna.
1Czyli: f0(t) = limh→0
f (t+h)−f (t)
8.3. ZBIEŻNOŚĆ METODY EULERA 63
Plan dowodu:
1. Szacujemy błąd metody – czyli odległość rozwiązania prawdziwego od przybliżo-nego.
2. Szacowanie powyższe prowadzi do nierówności rekurencyjnej. 3. Korzystamy z modułu ciągłości2, aby uprościć szacowanie.
4. Rozwiązujemy nierówność rekurencyjną, tzn. chcemy otrzymać szacowanie błędu metody w węźle tj, które nie będzie zależne od dobru j a jednie od h. Korzystamy z warunku Lipschitza, który spełnia funkcja f zgodnie z założeniami naszego twierdzenia.
5. Otrzymujemy oszacowanie, które pozwala w prosty sposób pokazać tezę, tzn. pokazujemy, że jeśli długość przedziału h dąży do zera, to błąd metody w j-tym węźle również.
Dowód. Wprowadźmy na początek oznaczenie. Niech ∆j := |yj− x(tj)| oznacza błąd
metody Eulera w j-tym węźle. Wówczas mamy:
∆j+1= |yj+1− x(tj+1)|(def ) = yj+ hf (tj, yj) − x(tj) − Z tj+1 tj f (s, x(s))ds = Skorzystaliśmy właśnie z definicji ∆j, yj oraz zapisaliśmy x(tj+1) w postaci całki (x jest rozwiązaniem zagadnienia, a nie przybliżeniem). Teraz do naszego wyrażanie pod wartością bezwzględną dodamy i odejmiemy hf (tj, x(tj)), aby potem zastosować pewne szacowanie. Mamy więc:
= yj− x(tj) + hf (tj, yj) − hf (tj, x(tj)) + hf (tj, x(tj)) − Z tj+1 tj f (s, x(s))ds
Zauważmy teraz, że h = tj+1− tj, a wyrażenie f (tj, x(tj)) oczywiście jest stałą (przy ustalonym j), tak więc hf (tj, x(tj)) = (tj+1 − tj)f (tj, x(tj)) = Rtj+1
tj f (tj, x(tj))ds. Ko-rzystając z tego, wyrażenie hf (tj, x(tj)) przeniesiemy pod całkę. Skorzystamy również z nierówności trójkąta: ¬ |yj− x(tj)| + h|f (tj, yj) − f (tj, x(tj))| + Z tj+1 tj f (tj, x(tj)) − f (s, x(s))ds ¬
Skorzystamy teraz z definicji ∆j, oraz z warunku Lipschitza aby uprościć zapis3:
¬ ∆j+ hL|yj − x(tj)| + Z tj+1 tj |f (tj, x(tj)) − f (s, x(s))| ds = = ∆j+ hL∆j + Z tj+1 tj |f (tj, x(tj)) − f (s, x(s))| ds =
2Pojęcie to zostanie wyjaśnione w ramach dowodu.
3Każdy z trzech składników naszej sumy zapisujemy w innej postaci: |yj− x(tj)| = ∆j, h|f (tj, yj) −
f (tj, x(tj))| ¬ hL|yj − x(tj)|, a w całce wchodzimy z wartością bezwzględną pod całkę (po takiej operacji wartość albo wzrośnie albo nie zmieni się).
64 ROZDZIAŁ 8. PRZYBLIŻANIE ROZWIĄZANIA – METODA EULERA
= (hL + 1)∆j +
Z tj+1
tj
|f (tj, x(tj)) − f (s, x(s))| ds. Otrzymaliśmy więc następujące szacowanie (∗):
∆j+1 ¬ (hL + 1)∆j+
Z tj+1
tj
|f (tj, x(tj)) − f (s, x(s))| ds. Rozważmy teraz funkcję:
g(s) = f (s, x(s)), g : [a, b] → R.
Jest to funkcja ciągła, a skoro jest określona na przedziale domkniętym, to jest jedno-stajnie ciągła, czyli:
∀>0∃δ>0∀t,t0∈[a,b]|t − t0| < δ ⇒ |g(t) − g(t0)| < .
Z powyższej definicji wynika, iż wybór δ uzależniony jest od wyboru . Można wobec tego mówić o funkcji (w pewnym sensie). Można też (czego nie uzasadnimy) mówić o sytuacji odwrotnej, tzn. o doborze dla dowolnej δ, czyli:
∀t,t0∈[a,b]|t − t0| < δ ⇒ |g(t) − g(t0)| < ω(δ),
gdzie ω : [0, b − a] → R+ oraz limδ→0+ω(δ) = 0. Taka funkcja ω nosi nazwę moduł
ciągłości funkcji g. Można ją skonstruować, jednak nie będziemy tego robić. Zakładamy po prostu, że ona istnieje.
Skorzystamy z modułu ciągłości w naszej nierówności. Nierówność (∗) można bo-wiem zapisać prościej:
∆j+1 ¬ ∆j(1 + hL) +
Z tj+1
tj
ω(h)ds.
Co można zapisać jeszcze prościej:
∆j+1 ¬ ∆j(1 + hL) + hω(h)ds. (∗∗)
Jest to ostateczna forma tej nierówności. Jak widać, jest to nierówność rekurencyjna, bo ∆j+1 zależy do ∆j. Postaramy się teraz „rozwikłać” tą nierówność, tzn. podać oszacowanie, które nie będzie rekurencyjne. Dla ustalenia uwagi obniżamy na początek indeksy o 1. Mamy:
∆j ¬ ∆j−1(1 + hL) + hω(h) ¬
Dla wyrażenia ∆j−1 rekurencyjnie możemy zastosować oszacowanie (∗∗):
¬ (1 + hL)((1 + hL)∆j−2+ hω(h)) + hω(h) = = (1 + hL)2∆j−2+ (1 + hL)hω(h) + hω(h) ¬ Podobnie jak poprzednio, szacujemy wyrażenie ∆j−2:
¬ (1 + hL)2((1 + hL)∆j−3+ hω(h)) + (1 + hL)hω(h) + hω(h) = = (1 + hL)3∆j−3+ (1 + hL)2hω(h) + (1 + hL)hω(h) + hω(h) ¬ ¬ . . . ¬ (1 + hL)j∆0+ (1 + hL)j−1hω(h) + . . . + (1 + hL)hω(h) + hω(h) =
8.3. ZBIEŻNOŚĆ METODY EULERA 65 = (1 + hL)j∆0+ hω(h) j−1 X i=0 (1 + hL)i Otrzymaliśmy więc: ∆j ¬ (1 + hL)j ∆0+ hω(h) j−1 X i=0 (1 + hL)i.
Oczywiście formalnie rzecz biorąc, potrzebny jest teraz dowód indukcyjny poprawności tego oszacowania. Łatwiej jednak udowodnić podany niżej lemat.
Lemat 8.3.3. Jeśli zachodzi ∆j+1 ¬ A∆j + B, to zachodzi również ∆j ¬ Aj∆0 +
BPj−1 i=0Ai.
Dowód lematu: Prosta indukcja po j. Krok I: dla j = 1 mamy ∆1 ¬ A∆0+ B, czyli teza i założenie mają tą samą postać. Więc lemat jest spełniony dla j = 1. Krok II: zakładamy że dla j lemat jest poprawny i pokażemy, że jest poprawny dla j + 1:
∆j+1 ¬ A(∆j) + B ¬ A(Aj∆0+ B j−1 X i=0 Ai) + B = = Aj+1∆0+ B j X i=1 Aj+ B = Aj+1∆0+ B j X i=0 Aj.
Czyli lemat jest również poprawny dla j + 1. Zgodnie z zasadą indukcji matematycznej, lemat jest prawdziwy.
Wróćmy do naszego oszacowania. Zauważmy, że wzór można jeszcze uprościć, ze względu na fakt iż ∆0 = 0. Dzieje się tak, ze względu na to, że w metodzie Eulera węzeł o numerze 0 dany jest w warunku początkowym zagadnienia, czyli nie jest on w żaden sposób przybliżony. Mamy więc:
∆j ¬ hω(h)
j−1
X
i=0
(1 + hL)i =
Stosując wzór na skończoną sumę ciągu geometrycznego mamy: = hω(h)1 − (1 + hL) j 1 − (1 + hL) = hω(h) 1 − (1 + hL)j −hL = = ω(h)1 − (1 + hL) j −L = ω(h) (1 + hL)j − 1 L . Mamy więc: ∆j ¬ (1 + hL) j − 1 L ω(h).
Skorzystamy teraz z nierówności:
∀x∈R(1 + x) ¬ ex.
Możemy więc napisać, że 1 + hL ¬ ehL. Mamy więc:
∆j ¬ e
hLj − 1
66 ROZDZIAŁ 8. PRZYBLIŻANIE ROZWIĄZANIA – METODA EULERA
Zauważmy, że ze względu na to, że ex jest rosnąca, to wyrażenie z prawej strony nie-równości osiąga największą wartość dla j = n. Mamy więc:
∆j ¬ e
hLn− 1
L ω(h).
Ze względu na to, iż h = b−an możemy napisać:
∆j ¬ e
(b−a)L
L ω(h).
Widać teraz, że jeśli przejdziemy do granicy h → 0+, to wyrażenie po prawej stronie nierówności dąży do zera (ze względu na funkcję ω(h)). Zgodnie z twierdzeniem o trzech ciągach, ∆j musi więc też zbiegać do zera, bo z jednej strony ograniczone jest przez zero, a z drugiej (zgodnie z naszym oszacowaniem) przez coś co zbiega do zera. Czyli rzeczywiście, jeśli liczba punktów podziału naszego odcinka rośnie, to h maleje, i błąd w każdym z punktów maleje. Metoda Eulera jest więc zbieżna.
Rozdział 9
Równania cząstkowe – równanie
ciepła
9.1 Wprowadzenie
W rozdziale zajmiemy się nieco równaniami cząstkowymi. Zauważmy, że dotychczas wszystko co robiliśmy z równaniami różniczkowymi dotyczyło równań, w których nie-wiadoma była funkcją jednej zmiennej (jedno, bądź wielo wymiarową). W równaniach o których będziemy mówić tutaj, niewiadoma, jest funkcją wielu zmiennych i pojawiają się pochodne cząstkowe tychże funkcji. Na początku podamy kilka przykładów równań cząstkowych, co pozwoli nam zapoznać się z podstawowymi oznaczeniami i terminami.
9.1.1 Przykłady równań cząstkowych
Definicja 9.1.1 (laplasjan). Laplasjanem funkcji u(x1, x2, . . . , xn) nazywamy:
∆u(x1, x2, . . . , xn) = ∂ 2u ∂x2 1 + ∂ 2u ∂x2 2 + . . .∂ 2u ∂x2 n .
Uwaga 9.1.2. Często mamy do czynienia z sytuacją, gdy pierwsza zmienna funkcji,
jest w jakiś sposób wyróżniona – na przykład oznacza czas. Wygodnie jest wówczas zdefiniować specjalną „wersje” laplasjanu dla funkcji u(t, x1, . . . , xn):
∆xu(t, x1, x2, . . . , xn) = ∂ 2u ∂x2 1 +∂ 2u ∂x2 2 + . . . ∂ 2u ∂x2 n
W wielu przypadkach, taki właśnie laplasjan będziemy oznaczać przez ∆u.
Przykład 9.1.3 (równanie Laplace’a). Równanie cząstkowe może mieć postać:
∆u = 0,
dla x ∈ Ω, z warunkiem: u = ϕ na ∂Ω co oznacza, że funkcja u przyjmuje wartość ϕ na ∂Ω czyli na brzegu zbioru Ω.
Przykład 9.1.4 (równanie drgania). Innym przykładem równania cząstkowego, może
być, stosowane w fizyce, równanie drgania:
∂2u
∂t2 − ∆xu = f.
68 ROZDZIAŁ 9. RÓWNANIA CZĄSTKOWE – RÓWNANIE CIEPŁA
Przykład 9.1.5 (równanie dyfuzji). Równanie dyfuzji, ma postać:
∂u
∂t − ∆xu = g.
W równaniu tym funkcja u oznacza temperaturę w punkcie x = (x1, . . . , xn) przestrzeni, w czasie t.
Przykład 9.1.6 (równanie ciepła). Równanie ciepła, którym będziemy się zajmować
dalej w tym rozdziale, to równanie postaci:
∂u
∂t − ∆xu = 0.