1 Wykład no 11

2 Metody numeryczne rozwiązywania równań różniczkowych cząstkowychMetoda różnic skończonych (siatek) Uwagi ogólne Dane równanie różniczkowe cząstkowe opisane operatorem L: w obszarze  i warunki brzegowe:

3 W metodach różnicowych poszukuje się tablicy wartości węzeł pomocniczy h=(hx,hy) hy Xk hx Parametr h charakteryzuje siatkę h x węzeł podstawowy W metodach różnicowych poszukuje się tablicy wartości przybliżonych uh rozwiązania dokładnego u na zbiorze izolowanych punktów Xk (k=1,2,...,Nh ) zwanym siatką. Punkty Xk są nazywane węzłami siatki. Równania służące do wyznaczania wartości przybliżonych nazywamy równaniami różnicowymi.

4 Dla równania różniczkowego:w obszarze  z warunkami brzegowymi: otrzymujemy jego odpowiednik różnicowy: Zakładając, że problem opisany równaniem różniczkowym ma jednoznaczne rozwiązanie, to równania różnicowe będą jego odpowiednikiem jeżeli są spełnione następujące warunki:

5 1. Układ równań różnicowych posiada jednoznaczne rozwiązanie:dla każdego dopuszczalnego h. 2. Zbieżność do rozwiązania dokładnego u. Oznacza to, że rozwiązanie uh powinno przy h 0 dążyć do rozwiązania dokładnego u. Dla określenia zbieżności jest koniecznym wprowadzenie odpowiednich przestrzeni funkcyjnych i norm w nich.

6 Wprowadzamy przestrzeń funkcyjną U z normą || ||U , do której należy rozwiązanie dokładne u. oraz przestrzeń Nh- wymiarową Uh z normą || ||Uh , której elementami są układy Nh liczb i do której należy rozwiązanie uh. Normy || ||U i || ||Uh winny być zgodne, tzn. ponieważ funkcja u(x) jest określona w węzłach podstawowych Xk siatki, to mówimy, że normy są zgodne jeżeli zachodzi: dla każdego uU.

7 Przykłady norm zgodnych:- zbiór węzłów wewnętrznych.

8 Wielkość: nazywamy błędem rozwiązania przybliżonego uh. Uh dąży do rozwiązania dokładnego u(x), jeżeli Jeżeli można znaleźć taką funkcję (h), że to mówimy, że zostało znalezione oszacowanie błędu.

9 3. Stabilność Różnicowe zagadnienie brzegowe jest stabilne, jeżeli istnieją takie liczby >0 i h0>0, że dla dowolnych h że ||fh||Fh< zadanie brzegowe: posiada jedno jednoznaczne rozwiązanie i zachodzi: gdzie C stała niezależna od h.

10 Twierdzenie wiążące stabilność i zbieżność:Jeżeli zadanie różnicowe jest przybliżeniem różniczkowego problemu brzegowego: i rozwiązanie uh jest stabilne wtedy zachodzi i rząd zbieżności w funkcji h jest taki sam jak rząd aproksymacji.

11 Zastępowanie pochodnych ilorazami różnicowymi na siatce prostokątnejhk k-1 i-1 i i+1

12 Druga pochodna dodająck+1 Druga pochodna dodając stronami: k hi hk k-1 i-1 i i+1

13 k+1 Dla pochodnej mieszanej k hi hk k-1 i-1 i i+1

14 Konstrukcja warunków brzegowych na siatcer(i,k,) i 1. Przeniesienie wartości: Przyjmujemy: lub

15 2. Interpolacja liniowa dla warunku brzegowego Dirichleta i,k i+1,k i-1,k Zakładając liniowy rozkład rozwiązania między sąsiednimi węzłami mamy: i dla x=xi+ mamy:

16 3. Interpolacja liniowa dla warunku brzegowego Neumanna.hi i-1,k i,k hk n i,k-1

17 Dla uproszczenia rozważań będą analizowane przypadkiRównania eliptyczne Dla uproszczenia rozważań będą analizowane przypadki dwuwymiarowe x=(x,y) Warunki brzegowe:

18 10 y(k) 9 8 7 p 6 hy j-1 j j+1 5 l 4 3 2 (0,0) hx 1 1 2 3 4 5 6 7 8 9 x(i)

19 Dla węzłów wewnętrznych będzie:p j-1 j j+1 l lub w formie macierzowej:

20 p hy j-1 j j+1 l (0,0) hx

21 styczna m p-1 p (xkp,ykp) normalna Przyjmując: otrzymujemy:

22 Uwaga dotycząca błędu obliczeń. Generalnie jeżeli węzły nie leżą na krzywej brzegowej i liczymy metodą przeniesienia wartości, to dokładność obliczeń jest rzędu h. Jeżeli węzły na krzywej brzegowej bądź wyliczamy wartości funkcji brzegowej interpolując liniowo dokładność wzrasta do h2. W zagadnieniu Dirichleta oprócz trudności z wyznaczeniem wartości brzegowych nie ma innych problemów i otrzymany układ równań algebraicznych najczęściej można rozwiązać bez kłopotów. Sytuacja może się komplikować przy zagadnieniu Neumanna. Siatki praktycznie nie stosowane w zagadnieniach eliptycznych liniowych i nieliniowych.

23 Równania paraboliczneRównania opisujące ewolucję układu w czasie Równania paraboliczne Dane jednowymiarowe równanie przewodnictwa: t k+1 k k-1  x h i-1 i i+1 1

24 t k+1 k k-1  x h i-1 i i+1 N Oznaczamy operator różnicy II rzędu: i wprowadzamy schematy różnicowe z wagą : i

25 Problem brzegowy jest aproksymowany przezdla i=1,2,...N dla k=1,2,...,K. Warunki zgodności

26 Schemat sześciowęzłowyk+1 k i i i+1 k+1 k Jeżeli =0 schemat jest nazywany jawnym lub explicite i i i+1 k+1 k Jeżeli 0 schemat jest nazywany niejawnym lub implicite i i i+1

27 Wartości w warstwie k+1 otrzymujemy rozwiązując układ równań: k+1 k i i i+1 k+1 k Czysto niejawny schemat: i k+1 k Schemat Cranka - Nicholsona: i i i+1

28 Oszacowanie dokładności aproksymacji.Rozwiązanie dokładne zagadnienia brzegowego jest u(x,t) i jego wartość w węzłach (xi,tk) siatki będzie oznaczana u(i,k). Rozwiązanie zagadnienia brzegowego sformułowanego dyskretnie jest Błąd aproksymacji jest

29 Dla oceny błędu w kroku k-tym wprowadza się normę, np.: lub wynika, że Z i podstawiając do w miejsce otrzymujemy równoważne zadanie różnicowe dla

30 gdzie błąd schematu różnicowego w stosunku do rozwiązania dokładnego u(x,t).

31 Mówimy, że przybliża rozwiązanie problemu brzegowego z dokładnością rzędu (m,n) lub O(hm+n), jeżeli spełnia nierówność: gdzie M - stała.

32 Dla uproszczenia zapisu wprowadzamy:i mamy: Uwzględniając powyższe równości i podstawiając do

33 Rozwijając funkcje u(n,p) w szereg Taylora w otoczeniu mamy Rozwijając funkcje u(n,p) w szereg Taylora w otoczeniu punktu: xi,tk+0.5  oraz wprowadzając oznaczenie: Będziemy mieli:

34 Uwzględniając powyższe zależności możnazapisać w postaci:

35 Ale gdyż u jest rozwiązaniem dokładnym i w każdym punkcie obszaru spełnia równanie paraboliczne. Uwzględniając, że mamy

36 Jeżeli wyrażenie w nawiasie kwadratowym jest równe zeru, tzn.:to schemat ma dokładność O(h4+2) oraz W obliczeniach numerycznych wygodniej przyjąć:

37 jest schematem o podwyższonej dokładności wynoszącej:Jeżeli  =0.5 jak w schemacie Cranka-Nicholsona, to

38 i dla =0.5 mamy: Dla zachowania oceny zbieżności O(h2+2) należy przyjąć: lub Jeżeli 0.5 i *, to dokładność obliczeń jest rzędu O(h2+).

39 Wesołych Świąt Wesołych Świat

40 Stabilność Zbadamy zachowanie się schematu: jawnego tj.  =0

41 Schemat jest Przyjmijmy: i załóżmy, że warunek początkowy w punkcie i-tym jest dany z błędem  . Badamy jak przenosi się błąd na siatce.

42 t x N Schemat jawny z =0 i

43 t x N Schemat jawny z =0 i Obliczenia z dokładnością do 2 miejsc

44 Analiza stabilności metodą spektralnąNiech rozwiązanie jednorodnego zagadnienia różnicowego: ma postać gdzie

45 ale Po podstawieniu do mamy

46 Po podzieleniu przez otrzymujemy równanie: Warunek konieczny stabilności Neumanna stwierdza, że schemat różnicowy jest stabilny, jeżeli Dla =0 mamy

47 i na podstawie kryterium Neumanna otrzymujemy:czyli Dla = znajdujemy warunek na stosunek: Warunkiem zbieżności schematu jawnego jest spełnienie powyższego warunku. Warunek jest również prawdziwy dla schematu jawnego w przypadku wielowymiarowego równania parabolicznego.

48 Dowolne >0. Ocenę prowadzimy przy =. czyli Prawa nierówność jest spełniona dla dowolnych , a z lewej mamy:

49 Dla  spełniających nierówność:warunek: jest spełniony dla dowolnego stosunku W szczególności schemat Cranka-Nicholsona =0.5 jest stabilny dla dowolnego stosunku kroków

50 Dla schematu o podwyższonej dokładnościmamy bo Przedstawione rozważania można rozszerzyć na przypadki wielowymiarowe jak również na równania o zmiennych współczynnikach. W przypadku równań wielowymiarowych ocena zbieżności zależy również od sposobu aproksymacji warunków brzegowych podobnie jak w przypadku równań eliptycznych.

51 Równania hiperboliczneJako przykład zostanie rozpatrzone równanie linii długiej bez strat o długości L. Dla napięcia u mamy: Wprowadzamy siatkę prostokątną:

52 t k+1 k k-1  x h i-1 i i+1 N i funkcję węzłową oznaczamy: Przyjmujemy aproksymację pochodnych: i

53 i rozpatrujemy następujący schemat trójwarstwowyz parametrem >0: gdzie warunek początkowy konstruujemy tak, aby zachować rząd aproksymacji O(2).

54 Mamy czyli Z równania falowego mamy: Warunek początkowy dla pierwszej pochodnej będzie określony z dokładnością O(h2+2), jeżeli przyjąć, że czyli

55 Ostatecznie schemat różnicowy dla rozwiązania równaniafalowego jest Ocena dokładności aproksymacji Postępujemy podobnie jak poprzednio, a więc niech

56 jest rozwiązaniem różnicowego zagadnieniagdzie a u(xi,tk) jest rozwiązaniem problemu brzegowego: w punkcie xi,tk.

57 Pisząc otrzymujemy: gdzie Z konstrukcji warunku początkowego dla pochodnej wynika, że

58 Rozwijając w szereg Taylora mamy:Korzystając z otrzymanego wyniku mamy: Stąd otrzymujemy z dokładnością do małych 4-go rzędu:

59 Podobnie z dokładnością do małych 4-go rzędu. Ostatecznie otrzymujemy ocenę błędu: Funkcja u spełnia równanie falowe, a więc stąd niezależnie od  !!!

60 Oznacza to również zależność odwrotną, a mianowicie dobór  zależy tylko i wyłącznie od stabilności, a nie ma wpływu na dokładność obliczeń. Stabilność Analizujemy stabilność schematu różnicowego przy jednorodnych warunkach brzegowych:

61 Przyjmujemy rozwiązanie w postaci:i podstawiając do równania różnicowego: po wykonaniu kilku przekształceń otrzymujemy: Dzieląc równanie przez k-1 i grupując wyrazy otrzymujemy równanie określające :

62 Badamy pierwiastki równania:gdzie przy =. Wyróżnik: Jeżeli <0, to równanie ma dwa sprzężone pierwiastki zespolone 1, 2 o module równym 1 co wynika ze wzoru Viety: Dla =0 otrzymujemy warunek Couranta:

63 który oznacza, że prędkość wędrówki fali na siatce h/ jest większa od prędkości fazowej. Przypadek 0 prowadzi do pierwiastków większych co do modułu od jedności i dlatego należy te przypadki odrzucić. Analizę można rozszerzyć na przypadki wielowymiarowe.

64 Wady i zalety metody różnicowej1. Proste konstruowanie siatki podziałowej. 2. Prosta konstrukcja układu równań różnicowych szczególnie w środowiskach izotropowych. 3. Opracowane oceny błędów metody i warunki stabilności. Wady: 1. Duże trudności z dobrą aproksymacją brzegu lub wymuszony mały krok siatki. 2. Trudności z utrzymaniem rzędu aproksymacji przy interpolacji warunków brzegowych. 3. Praktycznie konieczność obliczania całego obszaru z tym samym krokiem podziałowym.