Dyskretyzacja rwnania dyfuzji cd jawny Euler niejawny Euler

  • Slides: 70
Download presentation
Dyskretyzacja równania dyfuzji cd. jawny Euler niejawny Euler t +O(Dt) (błąd dyskretyzacji) t+Dt +O(Dt)

Dyskretyzacja równania dyfuzji cd. jawny Euler niejawny Euler t +O(Dt) (błąd dyskretyzacji) t+Dt +O(Dt) schemat Cranka-Nicolsona: +O(Dt 2) CN to odpowiednik wzoru trapezów dla dy/dt=f(t)

Dyskretyzacja równania dyfuzji cd. Schemat CN – niejawny, bez ograniczenia na krok czasowy ze

Dyskretyzacja równania dyfuzji cd. Schemat CN – niejawny, bez ograniczenia na krok czasowy ze względu na stabilność, jawny Euler niejawny Euler drugi rząd dokładności dla równania adwekcji – poznaliśmy schemat jawny tego samego rzędu dokładności, t t+Dt który powstaje przy inaczej wprowadzonej dyskretyzacji czasu. +O(Dt) (błąd dyskretyzacji) Schemat ten (leapfrog) był odpowiednikiem metody punktu pośredniego: +O(Dt) schemat Cranka-Nicolsona: +O(Dt 2) CN to odpowiednik wzoru trapezów dla dy/dt=f(t)

leap – frog (jawny, dwustopniowy) całkiem nieźle sprawdzał się dla równania adwekcji (równie dobrze

leap – frog (jawny, dwustopniowy) całkiem nieźle sprawdzał się dla równania adwekcji (równie dobrze co CN) czy zadziała dla równania dyfuzji? +O(Dx 2 )+O(Dt 2 ) dokładność jak CN 2 r analiza von Neumanna szukamy rozwiązań postaci: -

leap – frog analiza v. N cd mk =(1 -cos(k. Dx)) 0 dla bezwzględnej

leap – frog analiza v. N cd mk =(1 -cos(k. Dx)) 0 dla bezwzględnej stabilności |gk | ma być mniejsze od 1

leap – frog analiza v. N cd mk =(1 -cos(k. Dx)) 0 załóżmy, że

leap – frog analiza v. N cd mk =(1 -cos(k. Dx)) 0 załóżmy, że r małe rozwijamy g w szereg Taylora : rozwiązanie ogólne: właściwe równania dyfuzji pasożytnicze: rośnie co do modułu z n znak oscyluje z iteracji na iteracje rozwiązanie pasożytniczne pojawi się i doprowadzi do eksplozji leapfrog nie jest dobrym (bzwz. stabilnym) schematem dla r. dyfuzji

leapfrog: 2 r dla r=1/2 widzimy, że schemat jest symetryczny względem czasu licząc równanie

leapfrog: 2 r dla r=1/2 widzimy, że schemat jest symetryczny względem czasu licząc równanie wstecz dostaniemy ten sam przepis ale rozwiązanie równania dyfuzji NIE jest symetryczne względem czasu (t: =-t) - w przeciwieństwie do rozwiązania równania adwekcji W problemie stygnącego pręta: zanik temperatury wyznacza kierunek upływu czasu

odwrotny problem przewodnictwa cieplnego problem prosty : zadajemy warunki brzegowe oraz początkowe pytanie: co

odwrotny problem przewodnictwa cieplnego problem prosty : zadajemy warunki brzegowe oraz początkowe pytanie: co stanie się w przyszłości (tak wprowadzane są problemy w teorii równań różniczkowych) W praktyce, często chcemy znaleźć rozwiązania dla problemu odwrotnego: znamy obecny rozkład temperatury. Jaki był rozkład w przeszłości? Jakie były warunki brzegowe? Jaki był warunek początkowy ? [typowy problem pomiarów, nie tylko zależnych od czasu] warunki brzegowe u(x=0, t)=u(x=1, t)=0 problem: dane u(x, t=T) szukane: u(x, t=0)

O czasie i problemie odwrotnym. . . N=100, dx=1. 0/(N), D=1 dt=dx**2/d/2/10 (malutki) CN

O czasie i problemie odwrotnym. . . N=100, dx=1. 0/(N), D=1 dt=dx**2/d/2/10 (malutki) CN problem: T(x, t) = 1 wewnątrz T(x, t) = 0 na zewnątrz chcemy wrócić do warunku początkowego ustawiamy dt: =-dt x t liczymy do przodu potem chcemy wrócić

O czasie i problemie odwrotnym. . . N=100, dx=1. 0/(N), D=1 dt=dx**2/d/2/10 (malutki) CN

O czasie i problemie odwrotnym. . . N=100, dx=1. 0/(N), D=1 dt=dx**2/d/2/10 (malutki) CN problem: T(x, t) = 1 wewnątrz T(x, t) = 0 na zewnątrz chcemy wrócić do warunku początkowego ustawiamy dt: =-dt nieco dłużej = eksplozja x t liczymy do przodu potem chcemy wrócić i bum!

problem odwrotny do równania dyfuzji: wszystkie metody różnic skończonych okazują się niestabilne dla ujemnego

problem odwrotny do równania dyfuzji: wszystkie metody różnic skończonych okazują się niestabilne dla ujemnego kroku czasowego [r =DDt/Dx 2 < 0 ] wyprowadzony wcześniej z analizy von Neumanna warunek: stabilności bezwzględnej: (q=1/2 odpowiada CN) dla Dt <0 [r<0] warunek prawy nie jest spełniony schemat CN nie jest stabilny dla równania dyfuzji rozwiązywanego wstecz

niezależnie od startu rozkład T po pewnym czasie będzie miał kształt sin(px) Układ zapomina

niezależnie od startu rozkład T po pewnym czasie będzie miał kształt sin(px) Układ zapomina o warunku początkowym problem obiektywnie trudny odwracamy znak czasu: gdy tylko w wyniku niedokładności pojawi się składowa o wysokim n – natychmiast eksploduje Nie zawsze problem z cofaniem się wstecz w czasie jest trudny: dla równania adwekcji– jest równie łatwy jak początkowy (zmiana dt równoważna zmianie kierunku prędkości unoszenia v)

niezależnie od startu rozkład T po pewnym czasie będzie miał kształt sin(px) problem obiektywnie

niezależnie od startu rozkład T po pewnym czasie będzie miał kształt sin(px) problem obiektywnie trudny możliwe rozwiązanie: szukać warunków początkowych T(x, t=0), dla których jesteśmy najbliżej danych wejściowych [T(x, t=T)] rozwiązywać równanie dla dt>0 i porównywać wynik numeryczny dla t=T z zadanym rozkładem – co wymaga znacznie większego nakładu obliczeń niż w problem podstawowy

odwrotny problem przewodnictwa cieplnego opiszemy rozwiązanie warunków brzegowych u(x=0, t)=u(x=1, t)=0 problem: dane u(x,

odwrotny problem przewodnictwa cieplnego opiszemy rozwiązanie warunków brzegowych u(x=0, t)=u(x=1, t)=0 problem: dane u(x, t=T) szukane: u(x, t=0) policzone schematem CN dla N=100 dx=1. 0/(N) D=1 dt=dx**2/D/2 100 kroków czasowych Jeden z możliwych algorytmów – wykorzystuje liniowość równania

1) Jeden z możliwych algorytmów – wykorzystuje liniowość równania wybrać bazę niezależnych liniowo funkcji

1) Jeden z możliwych algorytmów – wykorzystuje liniowość równania wybrać bazę niezależnych liniowo funkcji g(x) określonych na przedziale (0, 1) np. gi (x) = (x-1/2)i 2) dla każdego warunku początkowego rozwiązać równanie przewodnictwa cieplnego do chwili T dostaniemy bazę funkcji hi(x) normalizujemy je tak aby (hi, hi)=1 zgodnie z tym warunkiem normalizujemy również gi

3) równanie jest liniowe ewolucja czasowa gi hi wyliczymy przybliżony warunek początkowy [wsp. d]

3) równanie jest liniowe ewolucja czasowa gi hi wyliczymy przybliżony warunek początkowy [wsp. d] jeśli rozłożymy rozwiązanie w chwili T w bazie funkcji hi rozłożyć: np. : metodą najmniejszych kwadratów

+ z niestety A bywa źle uwarunkowana bo hi mają tendencję do „upodabniania się”

+ z niestety A bywa źle uwarunkowana bo hi mają tendencję do „upodabniania się” nawet jeśli g bardzo różne niestety = raczej reguła dla problemów odwrotnych

Wyniki: dokładny warunek początkowy rozwiązanie problemu odwrotnego w bazie wielomianowej (i=0, 1, . .

Wyniki: dokładny warunek początkowy rozwiązanie problemu odwrotnego w bazie wielomianowej (i=0, 1, . . . 10) dokładny wynik: warunek początkowy był x(x-1)(x-1/4) baza dla i=0, 1, . . . 10

„upodabnianie się funkcji bazowych”- nie jesteśmy bez wpływu na uwarunkowanie problemu – możemy wybrać

„upodabnianie się funkcji bazowych”- nie jesteśmy bez wpływu na uwarunkowanie problemu – możemy wybrać bazę tak, aby efekt zminimalizować cos(npx) wielomiany gaussowska wielocentrowa g h t t t

Równanie adwekcji – dyfuzji (schematy jawne) występuje np. w mechanice płynów i pyłów w

Równanie adwekcji – dyfuzji (schematy jawne) występuje np. w mechanice płynów i pyłów w transporcie ciepła itd. D 0 Euler: przedni czasowy, centralne przestrzenne: schemat: bezwzględnie stabilny gdy czysta dyfuzja v=0 oraz r 1/2 : bezwzględnie niestabilny gdy czysta adwekcji D =0 : dla adwekcji widzieliśmy, że obecność niezerowego D stabilizuje schemat posortujmy wyrazy w powyższym równaniu względem indeksu siatki przestrzennej:

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli. .

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli. .

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli ½ r

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli ½ r |a | /2 aby schemat był stabilny: który efekt ma być dominujący: adwekcja czy dyfuzja ? ?

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli ½ r

równanie AD, schemat Eulera zgodnie z zasadą max: schemat będzie stabilny jeśli ½ r |a | /2 (przewaga dyfuzji) liczba Peclet’a (komórkowa liczba Reynoldsa) podobny wniosek otrzymamy dla normy euklidesowej stosując analizę von Neumanna 1) zauważmy – krok czasowy nie ma wpływu na stabilność jeśli prędkość unoszenia duża w porównaniu ze stałą dyfuzji: siatka przestrzenna będzie musiała być bardzo drobna. Zazwyczaj łatwiej zgodzić się na drobne dt niż na drobne dx ze względu na ograniczenia pamięciowe. 2) jeśli D=0 (czysta adwekcja) – schemat niestabilny

Dla równania adwekcji lepiej sprawdzał się schemat upwind : (dla a>0 ) [uwaga!, teraz

Dla równania adwekcji lepiej sprawdzał się schemat upwind : (dla a>0 ) [uwaga!, teraz v>0 wieje w prawo(inny znak v)] zasada max:

Dla równania adwekcji lepiej sprawdzał się schemat upwind [uwaga!, teraz v>0 wieje w prawo(inny

Dla równania adwekcji lepiej sprawdzał się schemat upwind [uwaga!, teraz v>0 wieje w prawo(inny znak v)] zasada max: r 0 (jest), r +a 0 (jest bo v>0) oraz 2 r+a 1 warunek znacznie mniej restrykcyjny niż dla Eulera bo: stabilność można zapewnić małym krokiem czasowym Dla dowolnej siatki ! Czy odnajdujemy znane warunki stabilności dla czystej dyfuzji i czystej adwekcji ?

Dla równania adwekcji lepiej sprawdzał się schemat upwind [uwaga!, teraz v>0 wieje w prawo(inny

Dla równania adwekcji lepiej sprawdzał się schemat upwind [uwaga!, teraz v>0 wieje w prawo(inny znak v)] zasada max: r 0 (jest), r +a 0 (jest bo v>0) oraz 2 r+a 1 warunek znacznie mniej restrykcyjny niż dla Eulera bo: stabilność można zapewnić małym krokiem czasowym Dla dowolnej siatki ! odnajdujemy znane warunki stabilności dla czystej dyfuzji i czystej adwekcji

problemy z przewagą adwekcji i v zmieniającym znak (a zależne od położenia) v>0 v<0

problemy z przewagą adwekcji i v zmieniającym znak (a zależne od położenia) v>0 v<0 co, można zapisać jednym wzorem: (z uniknięciem instrukcji warunkowej)

problemy z przewagą adwekcji i v zmieniającym znak (a zależne od położenia) v>0 v<0

problemy z przewagą adwekcji i v zmieniającym znak (a zależne od położenia) v>0 v<0 co, można zapisać jednym wzorem: z tzw. schemat z różniczkowaniem pod wiatr uwaga: w schemacie upwind: czynnik dyfuzji wzrasta o extra |a|/2 (pojawia się dyfuzja numeryczna) (w centralnym ilorazie sztucznej dyfuzji nie ma i to jak widzieliśmy powód niestabilności schematu dla czystej adwekcji) centralny (bez numerycznej dyfuzji) :

Przykład: problem z przewagą adwekcji D=0. 01, v=1 warunek początkowy: u=1/2 dla x<1/2 x

Przykład: problem z przewagą adwekcji D=0. 01, v=1 warunek początkowy: u=1/2 dla x<1/2 x t rozwiązanie dokładne dyfuzja: widoczna w lekkim zaokrągleniu nieciągłości dla t>0 upwind dt=0. 025, dx=0. 05 a=0. 5, r=0. 1 widać znacznie przesadzoną dyfuzję iloraz centralny (bezwzględnie niestabilny) widać generację niestabilności (antydyfuzja = zaostrzanie kantów) aby zniwelować dodatkową (numeryczną dyfuzję) dla schematu upwind - mniejszy krok czasowy czy mniejszy krok przestrzenny ?

Nieliniowe równania paraboliczne Dla równań liniowych (np. dyfuzji, dyfuzji+adwekcji) schematy jawne sprowadzają się do

Nieliniowe równania paraboliczne Dla równań liniowych (np. dyfuzji, dyfuzji+adwekcji) schematy jawne sprowadzają się do wykonania wielu podstawień w każdym kroku niejawne prowadzą do układu równań liniowych. Zastanowimy się jak rozwiązać równanie nieliniowe. schemat niejawny, jednopoziomowy, centralne przestrzenne przedni czasowy, ważona prawa strona (dla q=1/2 - CN), weźmy nieliniowe równanie dyfuzji na m=1 się już znamy

u(x, t=0)=exp(-x 2/25), pudło (-30, 30), Dx=1, Dt=. 1 dyfuzja nieliniowa m=5 m=1 zwykła

u(x, t=0)=exp(-x 2/25), pudło (-30, 30), Dx=1, Dt=. 1 dyfuzja nieliniowa m=5 m=1 zwykła dyfuzja CN

warunek początkowy oraz niejednorodność w chwili początkowej = do wyjaśnienia różnic w rozwiązaniu m=1

warunek początkowy oraz niejednorodność w chwili początkowej = do wyjaśnienia różnic w rozwiązaniu m=1 u m=5 Prawa strona równania mówi tutaj: „rośnij!” „nie zmieniaj się!” umxx „malej!” widzimy, że krańce pakietu = bez zmian. błyskawiczne stłumienie maksimum, wyrównanie brzegów

Nieliniowe równania paraboliczne zapiszemy jako układ równań nieliniowych dla q=0 – jawny schemat –

Nieliniowe równania paraboliczne zapiszemy jako układ równań nieliniowych dla q=0 – jawny schemat – nadal forma podstawieniowa (nawet dla nieliniowego równania) dla q 0 – schemat niejawny – metoda Newtona lub iteracja funkcjonalna

CN + iteracja funkcjonalna pierwszy krok czasowy, uzgodnienie punktu w centrum x=0 m=5 Dt=.

CN + iteracja funkcjonalna pierwszy krok czasowy, uzgodnienie punktu w centrum x=0 m=5 Dt=. 1 m=1 Dt=. 1

nieliniowe równanie dyfuzji CN, zbieżność iteracji funkcjonalnej punkt centralny, pierwszy krok czasowy m=5 Dt=.

nieliniowe równanie dyfuzji CN, zbieżność iteracji funkcjonalnej punkt centralny, pierwszy krok czasowy m=5 Dt=. 2 m=1 Dt=. 2 widzimy, że iteracja funkcjonalna nie rokuje dobrze dla zbieżności równania nieliniowego przy dłuższym kroku czasowym

jeśli z iteracją kłopoty może zastosować schemat jawny zamiast CN ? Dt=0. 3, 100

jeśli z iteracją kłopoty może zastosować schemat jawny zamiast CN ? Dt=0. 3, 100 kroków samouzgodnienia iteracją funkcjonalną CN: jawny: Zaczyna się dobrze pojawiają się wartości 1014 po czym pakiet zanika

A może schemat jawny zamiast CN ? Dt=0. 3, 100 iteracji CN schemat jawny

A może schemat jawny zamiast CN ? Dt=0. 3, 100 iteracji CN schemat jawny pojawiają się wartości 1014 po czym pakiet znika 1) niejawność schematu jest potrzebna 2) iteracja funkcjonalna się nie sprawdza metoda Newtona

metoda Newtona dla nieliniowego równania dyfuzji rozwiązania schematu dla n+1 kroku czasowego spełniają układ

metoda Newtona dla nieliniowego równania dyfuzji rozwiązania schematu dla n+1 kroku czasowego spełniają układ równań nieliniowych: F(Un+1)=0 rozwinięcie Taylora przybliżony wektor Un+1 w k-tej iteracji n+1 – znaczy n+1 chwila czasowa układ równań liniowych na poprawę przybliżenia Vk+1: =Vk+(Un+1 -Vk)

metoda Newtona dla nieliniowego równania dyfuzji U 0=UJ=0 macierz Jakobiego J-1 m. Jakobiego: trójprzekątniowa

metoda Newtona dla nieliniowego równania dyfuzji U 0=UJ=0 macierz Jakobiego J-1 m. Jakobiego: trójprzekątniowa

Wyniki [CN] dla pierwszego kroku czasowego m=1 Dt=. 1 iteracja funkcjonalna m=5 Dt=. 1

Wyniki [CN] dla pierwszego kroku czasowego m=1 Dt=. 1 iteracja funkcjonalna m=5 Dt=. 1 Metoda Newtona: 1 1 2 0. 970743147366556 2 0. 9922461168083 3 0. 970376491139719 3 0. 992246116808

Wyniki [CN] dla pierwszego kroku czasowego m=5 iteracja funkcjonalna Dt=. 2 Metoda Newtona: 1

Wyniki [CN] dla pierwszego kroku czasowego m=5 iteracja funkcjonalna Dt=. 2 Metoda Newtona: 1 1 m=1 Dt=. 2 Metoda Newtona: Równanie liniowe = zbieżność metody Newtona w jednej iteracji 2 0. 949526520122893 1 1 3 0. 947925874533601 2 . 98466247689 4 0. 947923849482469 3 . 98466247689

m=5, dt=1 z iteracją Newtona Wniosek: aby rozwiązać równania nieliniowe z rozsądnym krokiem czasowym

m=5, dt=1 z iteracją Newtona Wniosek: aby rozwiązać równania nieliniowe z rozsądnym krokiem czasowym potrzebna jest metoda niejawna do rozwiązania nieliniowych równań schematu - iteracja Newtona

Szacowanie błędów dla równań cząstkowych zależnych od czasu na przykładzie równania adwekcji czasowa i

Szacowanie błędów dla równań cząstkowych zależnych od czasu na przykładzie równania adwekcji czasowa i przestrzenna pochodna zastąpione przednim ilorazem różnicowym (jest to upwind dla v<0) z góry wiemy, że wyliczone wartości będą różnić się od wartości dokładnych o pewną wartość zależną od pochodnych rozwiązania dokładnego szacowanie błędów: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań metod o innym rzędzie dokładności -ale w praktyce ta wiedza nie przyda nam się do ilościowego oszacowania popełnionego błędu

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań metod o innym rzędzie dokładności błędy lokalne dwóch metod (zakładamy, że lokalny błąd czasowy jest o 1 wiekszy niz przestrzenny) np: dla p=1, pierwsza metoda: U – upwind [O(dt 2), O(dx)], druga: V - CN[O(dt 3), O(dx 2)]

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań metod o innym rzędzie dokładności błędy lokalne dwóch metod (zakładamy, że lokalny błąd czasowy jest o 1 wiekszy niz przestrzenny) np: dla p=1, pierwsza metoda: U – upwind [O(dt 2), O(dx)], druga: V - CN[O(dt 3), O(dx 2)] błąd dokładniejszego schematu: zaniedbywalny w porównaniu z błędem mniej dokładnego różnica V oraz U daje oszacowanie błędu gorszego schematu strategia: do ewolucji czasowej używamy V, możemy wypowiedzieć się o błędzie U

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań

szacowanie a posteriori: 2 opcje 1) porównanie rozwiązań na różnych siatkach 2) porównanie rozwiązań metod o ekstrapolacja Richardsona innym rzędzie dokładności używamy jednego schematu lecz dwóch siatek: (Dx, Dt) , oraz (Dx/2, Dt/2) t(n+1/2) t (n) t(n) x (j) Błąd w chwili n+1 w punktach rzadkiej siatki mamy:

mamy oszacowanie błędu obydwu rozwiązań ale tylko na rzadkiej siatce. . . co dla

mamy oszacowanie błędu obydwu rozwiązań ale tylko na rzadkiej siatce. . . co dla automatycznej kontroli Dt całkowicie wystarczy

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład upwind (iloraz przedni czasowy, wsteczny przestrzenny)

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład upwind (iloraz przedni czasowy, wsteczny przestrzenny) dokładność ilorazów przestrzennych i czasowych identyczna (p=1) błąd lokalny O(Dx)+O(Dt 2) v=1 warunek początkowy: u(x)=sin(x) rozwiązanie dokładne u(x, t)=sin(x-t) x 0 2 p zawsze u(0)=u(2 p) = zastosujemy periodyczne warunki brzegowe

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład xj=(j-1)Dx, Dx=2 p/J j=1, . .

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład xj=(j-1)Dx, Dx=2 p/J j=1, . . . J, u(x, t=0) J=4 Dx=p/2 Dt=p/4 u(x, t=Dt) po dwóch krokach Dt’=Dt/2 u U(J=8) J’=8 Dx’=Dx/2 Dt’=Dt/2 gdzie się pokrywają: szacujemy błąd po jednym Dt U(J=4) t szacujemy błąd -0. 1035533365 0. 1035533983 0. 1035533892 -0. 1035534073 błąd faktyczny -0. 1035534603 0. 1035533983 0. 1035533892 -0. 1035534073 rewelacyjny wynik

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład oszacowanie błędu w jednym kroku Dt

ekstrapolacja Richardsona dla równań różniczkowych cząstkowych - przykład oszacowanie błędu w jednym kroku Dt bardzo dokładne: wykorzystać do poprawy dokładności gęsta siatka: niebieska algorytm: w chwili tn znamy wartości funkcji na gęstej siatce 1) przepisujemy je naprzemiennie na dwie rzadsze siatki: czerwoną i czarną 2) wykonujemy krok Dt dla każdej z nich 3) wykonujemy dwa kroki Dt/2 na gęstszej siatce 4) szacujemy i odcinamy błędy w kroku t+Dt

upwind z ekstrapolacją Richardsona i usunięciem błędu dokładny t gęstsza: J=8 Dx’=Dx/2 Dt’=Dt/2 x

upwind z ekstrapolacją Richardsona i usunięciem błędu dokładny t gęstsza: J=8 Dx’=Dx/2 Dt’=Dt/2 x x upwind: bez obcięcia błędu = rozwiązanie zanika (dyfuzja)

Równanie dyfuzji oraz dyfuzji-adwekcji – typowe paraboliczne (opisuje dążenie do równowagi). Dziś zajmiemy się

Równanie dyfuzji oraz dyfuzji-adwekcji – typowe paraboliczne (opisuje dążenie do równowagi). Dziś zajmiemy się typowym równaniem hiperbolicznym (oscylacje: mechaniczne, elektryczne, elektro-magnetyczne) równanie falowe (dla struny) u(x, t) x=l x=p x+dx a b T 2 x T 1 siła naciągu struny T (kierunek poziomy): (II zasada Newtona F=ma)

Równanie dyfuzji oraz dyfuzji-adwekcji – typowe paraboliczne (dążenie do równowagi). Dziś zajmiemy się typowym

Równanie dyfuzji oraz dyfuzji-adwekcji – typowe paraboliczne (dążenie do równowagi). Dziś zajmiemy się typowym równaniem hiperbolicznym równanie falowe (dla struny) (oscylacje) u(x, t) x=l x=p x+dx a b uwaga: T 2 x T 1 dx 0 (prędkość rozchodzenia się drgań)

c – stałe: Ogólne rozwiązanie dla nieskończonego ośrodka (d’Alamberta) dowolna funkcja drgania rozchodzące się

c – stałe: Ogólne rozwiązanie dla nieskończonego ośrodka (d’Alamberta) dowolna funkcja drgania rozchodzące się bez zmiany kształtu [brak dyspersji w równaniu falowym] Liniowość równania: zasada superpozycji

Liniowość równania i zasada superpozycji: Sygnały rozchodzą się niezależnie od siebie F=exp(-(x-0. 5+ct)2) +exp(-(x+0.

Liniowość równania i zasada superpozycji: Sygnały rozchodzą się niezależnie od siebie F=exp(-(x-0. 5+ct)2) +exp(-(x+0. 5 -ct)2) Sygnały mijają się bez zmiany kształtu [(jedna fala przenika drugą. ] ponieważ równanie liniowe: jeśli wskażemy bazę zupełną funkcji ze znaną ewolucją czasową = problem rozwiązany baza: mody normalne (fale stojące) (drgania własne)

baza: mody normalne (fale stojące) (drgania własne) Dwupunktowe warunki brzegowe u(0, t)=u(L, t)=0 x=L

baza: mody normalne (fale stojące) (drgania własne) Dwupunktowe warunki brzegowe u(0, t)=u(L, t)=0 x=L x=0 Poszukajmy rozwiązań, w których tylko amplituda (a nie kształt fali) nie zależy od czasu: u(x, t)=X(x)T(t) t=0 t=t 2 t=t 3 t=t 1 T(t)=cos(w t+f)= C cos(w t)+D sin(w t) [gdy gęstość struny zmienna c może być funkcją położenia]

Równanie na część przestrzenną fal stojących (drania własne, drgania normalne) Dla c niezależnego od

Równanie na część przestrzenną fal stojących (drania własne, drgania normalne) Dla c niezależnego od x: k-liczba falowa, wektor falowy k = 2 p/ l tutaj l długość fali k=w / c Xn(x)=sin(knx) WB: spełnione, gdy X(0)=X(L)=0 kn=np /L Fale stojące: Między warunkami brzegowymi całkowita liczba połówek długości fal. . 0 L

warunki brzegowe: kwantyzacja k kwantyzacja w Xn(x)=sin(knx), Tn=sin(wnt) , cos(wnt) kn=np /L k=w /

warunki brzegowe: kwantyzacja k kwantyzacja w Xn(x)=sin(knx), Tn=sin(wnt) , cos(wnt) kn=np /L k=w / c T T oznacza naciąg struny wn =ckn przestrzenne drgania własne nie zależą od c, ale częstości tak. Wiemy, że niższe tony dają struny o większej grubości [ r ]. Wiemy również, że im silniej struna naciągnięta tym wyższy dźwięk.

Drgania własne dla zmiennej gęstości struny W przypadku ogólnym [c=c(x)] przyda się rachunek numeryczny.

Drgania własne dla zmiennej gęstości struny W przypadku ogólnym [c=c(x)] przyda się rachunek numeryczny. Wyliczyć Xn oraz wn r(x) Dyskretyzujemy drugą pochodną, liczymy Xn(x+dx)

Równanie własne z warunkami brzegowymi: Metoda strzałów. w – parametr równania w - dokładna

Równanie własne z warunkami brzegowymi: Metoda strzałów. w – parametr równania w - dokładna wartość własna w<w w=w w>w w X 0 wstawić warunek brzegowy ale co wstawić za X 1 ? ? (dla równania Poissona opisywanego dla metody Numerowa to był poważny problem) dla drgań własnych wstawiamy cokolwiek (równanie własne jest jednorodne rozwiązania określone co do stałej multiplikatywnej )

Test metody dla r(x)=1 (L=1, T 0=1) Analityczne: kn=np /L Miejsca zerowe – wartości

Test metody dla r(x)=1 (L=1, T 0=1) Analityczne: kn=np /L Miejsca zerowe – wartości własne przy których funkcje własne spełniają prawy warunek brzegowy X 0 x x 1 0 1 2 3 4 5 6 7 8 9 10 11 12

Przykład: r(x)=1+4 a(x-1/2)2 (struna cięższa przy mocowaniach) 4 3 2 1 0 0 a=0

Przykład: r(x)=1+4 a(x-1/2)2 (struna cięższa przy mocowaniach) 4 3 2 1 0 0 a=0 – częstości własne równoodległe Częstości własne maleją z a (cięższa struna) duże a – częstości grupują się w pary 1 W każdej parze: funkcja parzysta i nieparzysta. Środek struny – prawie nieważki, na częstości wpływ ma kształt funkcji przy brzegach – a tam zbliżony dla każdej funkcji z pary

Drgania własne a ogólne rozwiązania równania falowego Równanie ogólne: Warunki początkowe: u(x, t=0) oraz

Drgania własne a ogólne rozwiązania równania falowego Równanie ogólne: Warunki początkowe: u(x, t=0) oraz v(x, t=0)=du/dt Zadać wychylenie i prędkości rozłożyć warunki początkowe na drgania własne problem zależności czasowych jest rozwiązany

Drgania normalne a ogólne rozwiązania równania falowego Równanie ogólne: Warunki początkowe: u(x, t=0) oraz

Drgania normalne a ogólne rozwiązania równania falowego Równanie ogólne: Warunki początkowe: u(x, t=0) oraz v(x, t=0)=du/dt Zadać wychylenie i prędkości rozłożyć warunki początkowe na drgania własne problem zależności czasowych jest rozwiązany w chwili t=0, za kształt struny odpowiadają współczynniki cn a za prędkość – współczynniki sn

Superpozycja drgań własnych: Warunki początkowe Uwaga: dla równań dyfuzji i adwekcji warunek początkowy był

Superpozycja drgań własnych: Warunki początkowe Uwaga: dla równań dyfuzji i adwekcji warunek początkowy był tylko jeden czasowy rząd równania był = 1. Dla równania drugiego rzędu w czasie, wartość u(x, t=0) nie wystarczy dla jednoznacznego określenia rozwiązania. dla drgań własnych jednorodnej struny: Dyskretna sinusowa transformata Fouriera

rozkład na mody normalne na przedziale (0, L) Rozwinięcie w szereg Fouriera: g(x) =

rozkład na mody normalne na przedziale (0, L) Rozwinięcie w szereg Fouriera: g(x) = okresowa, odcinkowo ciągła z okresem T: Rozkład na drgania normalne a szereg Fouriera: drgania podległe warunkom brzegowym g(0)=g(L)=0. dla naszego problemu L to długość struny i nie ma interpretacji okresu (na strunie mieści się połowa długości fali).

Warunki Dirichleta zbieżności szeregu Fouriera Rozwinięcie Fouriera zbieżne w sensie jednorodnym N o ile

Warunki Dirichleta zbieżności szeregu Fouriera Rozwinięcie Fouriera zbieżne w sensie jednorodnym N o ile g(x) 1) całkowalna w kwadracie 2) odcinkowo ciągła rozwinięcie Fouriera dąży do g(x) „prawie wszędzie” tzn. poza punktami dyskretnymi punktami (rozwinięcie Fouriera jest wszędzie ciągłe!) Twierdzenie Dirichleta: W punktach nieciągłości szereg Fouriera zbieżny do g(x)=[ g(x-0)+g(x+0) ] / 2 tw. Dirichleta nie rozwiązuje wszystkich problemów

dla struny: pewien praktyczny problem z kanciastymi (nieróżniczkowalnymi) warunkami początkowymi. 1 Fala prostokątna -1

dla struny: pewien praktyczny problem z kanciastymi (nieróżniczkowalnymi) warunkami początkowymi. 1 Fala prostokątna -1 -p 0 p W punkcie nieciągłości = [g(0 -)+g(0+) ]/2 = (-1 + 1) / 2 =0

N=5 N=15 N=55 Nad nieciągłością wartość schodka przestrzelona o około 18% -p 0 p

N=5 N=15 N=55 Nad nieciągłością wartość schodka przestrzelona o około 18% -p 0 p 1+2 w Na PC pracujemy ze skończonymi bazami: N=15 N=55 N=100 Zjawisko Gibbsa w=0. 08949 (stała Wibrahama-Gibbsa) Równania różniczkowego przez rozkład warunku początkowego na drgania własne nie rozwiążemy dokładnie, jeśli ten jest nieciągły.

Na PC pracujemy ze skończonymi bazami. . . Zbieżność szeregu Fouriera w sensie bezwzględnym

Na PC pracujemy ze skończonymi bazami. . . Zbieżność szeregu Fouriera w sensie bezwzględnym Szereg jest bezwzględnie zbieżny jeśli można go obciąć na pewnym wyrazie rozwinięcia: Rozwinięcie fali prostokątnej nie jest bezwzględnie zbieżne: Bo ogólny szereg harmoniczny jest rozbieżny Wniosek: w skończonej bazie funkcji własnych możemy rozwiązywać tylko problemy z warunkiem początkowym, którego rozwinięcie w szereg Fouriera jest bezwzględnie zbieżne