Metody numeryczne - if.pw.edu.pl agatka/numeryczne/wyklad_09.pdf · Rozwa Ŝmy równanie struny zaczepionej…

Published on
28-Feb-2019
View
212
Download
0

Embed Size (px)

Transcript

Dr Piotr Fronczak Metody numeryczne Wykad nr 9 Rwnania rniczkowe zwyczajne - problemy brzegowe (BVP) Dotychczas omawialimy problemy pocztkowe rwania rniczkowe, w ktrych dane byy wartoci zmiennych zalenych (lub ich pochodne) dla pewnej szczeglnej wartoci zmiennej niezalenej. Teraz naszym zadaniem bdzie wyznaczenie spord funkcji speniajcych dane rwnanie rniczkowe zwyczajne, zdefiniowanych w rozwaanym obszarze, tych, ktre speniaj dodatkowe warunki na brzegu tego obszaru. Warunki takie nazywane s warunkami brzegowymi i s naoone na wartoci funkcji i jej pochodnych w wicej ni jednym punkcie tego obszaru. BVP s zwykle szczeglnym przypadkiem rwna rniczkowych czstkowych, ktrych rozwizaniem s funkcje czasu i pooenia, np. pole elektryczne, rozkad temperatury, prdko przepywu itp. 2 2 x y t y = rwnanie dyfuzji 2 2 2 2 x y t y = rwnanie falowe Jeli y = y(x) (nie zaley od czasu stany ustalone, rwnowagi), to otrzymujemy ogln posta BVP (drugiego rzdu): ),(2 2 dx dyyf dx yd=Dane jest rwnanie w dziedzinie bxa oraz okrelone s w pewien sposb warunki brzegowe. Typowe formy warunkw brzegowych Warunki brzegowe Dirichleta Dwie wartoci y(x) s dane jedna dla x = a, druga dla x = b. y(a) = Ya oraz y(b) = Yb Warunki brzegowe Neumanna Dwie wartoci dy/dx s dane jedna dla x = a, druga dla x = b. a ax Ddx dy= = b bx Ddx dy= = oraz Mieszane warunki brzegowe Robina a ax Caycdx dyc =+ = )(21 bbx Cbycdx dyc =+ = )(43oraz 4321 ,,, cccc - stae Metoda strzaw bxaxFyxkdx yd=+ );()( 2 2 Rozpatrzmy oglne rwnanie W celu jego zdyskretyzowania przyjmijmy Nabh /)( = gdzie N jest liczb punktw, na ktre dzielimy przedzia [a,b]. Dyskretyzujc drug pochodn )(2 2 2 11'' hOh yyyy iiii + += + iiiiii Fyk h yyy=+ + +2 11 2 otrzymujemy rwnanie iiiiii Fyk h yyy=+ + +2 11 2 Chcemy scakowa to rwnanie od x0 = a do xN = b, wic przedstawmy je w postaci iiiiii Fhykhyyy22 11 2 ++= + 1 2 11 2 102 2 Fhykhyyy ++=Czyli dla i = 1 Mamy dane y0 = y(a) = 0, ale y1 jest nieznane. Znajomo y1 jest rwnoznaczna ze znajomoci y dla x = 0: h yyy iii +1' 001 yhyy + Rozwamy rwnanie struny zaczepionej na obu kocach. Jedyn si dziaajc na element struny jest sia napre T. Jej warto w kierunku pionowym ( ))sin()sin()sin()sin( 1 iiTTTF =+= + Zakadajc, e kty s mae x yytg iiii = +++ 111 )()sin( oraz, e przypieszenie elementu struny jest proporcjonalne do wychylenia ya2= dostajemy += + x yyyTyxma iii 11 2 2)( 0)( 2 22 11 =+ + + yx yyyT iii Dla 0x 0'' 2 =+ yTy 0'' =+ yy T 2 = 2 2 2 2 x y t y = 2 2 dt yda = Zakadajc e otrzymalibymy 0)(,0)0(0'' ===+ Lyyyy Gdy > 0, to istnieje rozwizanie postaci )sin()cos( xBxAy += Uwzgldniajc warunki brzegowe otrzymujemy nieskoczenie wiele rozwiza ...,3,2,1,sin)( = = n L xnxyn Gdy = 0BAxy += Uwzgldniajc warunki brzegowe otrzymujemy rozwizanie trywialne y = 0. z wartociami wasnymi...,3,2,1, 2 22 == nL nn Poniewa problem fizyczny nie ma ujemnych wartoci wasnych, nie musimy analizowa przypadku < 0. Zatem mamy rozwizanie ...,3,2,1,sin)( = = n L xnxyn rwnania 0)(,0)0(0''2 22 ===+ LyyyL ny Wybierzmy n = 4 i L = 1. ( ) 0)1(,0)0(,4sin)( === yyxxy ( )xxy 4cos4)(' =Pochodna wynosi =0 =0 )162(4 222 hhy = iiiiii Fhykhyyy22 11 2 ++= + Oraz dalsze kroki zgodnie ze wzorem ze slajdu nr 6 )4cos(4 001 xhyy += 1 2 11 2 102 2 Fhykhyyy ++= A kolejne kroki rozwizania x yteor ynum 0 0.0000 0.0000 0.1 0.9511 1.2566 0.2 0.5878 0.5289 0.3 -0.5878 -1.0341 0.4 -0.9511 -0.9641 0.5 0.0000 0.6283 0.6 0.9511 1.2285 0.7 0.5878 -0.1113 0.8 -0.5878 -1.2753 0.9 -0.9511 -0.4255 1 0.0000 1.0963 Podzielmy obszar rozwiza [0,1] na 10 rwnych czci (h = 0.1) x yteor ynum 0 0.0000 0.0000 0.05 0.5878 0.6283 0.1 0.9511 1.0086 0.15 0.9511 0.9907 0.2 0.5878 0.5817 0.25 0.0000 -0.0570 0.3 -0.5878 -0.6731 0.35 -0.9511 -1.0235 0.4 -0.9511 -0.9699 0.45 -0.5878 -0.5333 0.5 0.0000 0.1138 0.55 0.5878 0.7160 0.6 0.9511 1.0355 0.65 0.9511 0.9462 0.7 0.5878 0.4834 0.75 0.0000 -0.1703 0.8 -0.5878 -0.7567 0.85 -0.9511 -1.0444 0.9 -0.9511 -0.9198 0.95 -0.5878 -0.4321 1 0.0000 0.2262 -1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000 0 0.2 0.4 0.6 0.8 1 1.2 yteor ynum -1.5000 -1.0000 -0.5000 0.0000 0.5000 1.0000 1.5000 0 0.2 0.4 0.6 0.8 1 1.2 yteor ynum Podzielmy obszar rozwiza [0,1] na 20 rwnych czci (h = 0.05) Zwykle nie znamy wartoci i wektorw wasnych ukadu. Trzeba je zgadn. Przeksztamy oglny problem brzegowy drugiego rzdu badx dyYbyYaybxayxf dx yd=== )(,)(),,( 2 2 zdla w ukad dwch rwna pierwszego rzdu: aYaywdx dy== )( pocz. war. z ),,( wyxfdx dw= BEZ WARUNKU POCZTKOWEGO Musimy znale warunek pocztkowy ===axdx dyaw )( Innymi sowy musimy znale nachylenie OKkrzywej y w punkcie a. Ya Yb x=a x=b2 < OK < 1 0)(,0)0(016'' 2 ===+ Lyyyy Wrmy do naszego przykadu (dla n=4, L=1) Przepiszmy to rwnanie w postaci ukadu dwch rwna pierwszego rzdu: == == )0(16 0)0( 2wyw' ywy' - parametr ukadu. Musimy znale miejsce zerowe funkcji bdu 0)1()1()( == yyE Zwykle powyszy ukad rwna bdziemy rozwizywa jedn z metod podanychna poprzednim wykadzie (np. Rungego-Kutty), ale na razie, korzystajc z metod analitycznych, zauwamy e Problem pocztkowy! 4 )4sin()( xCxy = A zatem nie znajdziemy staej C z warunku y(1) = 0. 0.2 0.4 0.6 0.8 1.0 -1.0 -0.5 0.5 1.0 Wybierzmy zatem na potrzeby dydaktyki inny warunek brzegowy: 587.0)7.0( =yMetoda bisekcji. 100 100 1 0 = = -100 -50 50 100 -6 -4 -2 2 4 0)(2 110 =+= r rrEE = 00 0)()( -10.5 -10.0 -9.5 -9.0 -0.06 -0.04 -0.02 0.02 0.04 2.10=r Sprawdmy: 17.10)7.04cos(4)'7.04sin( == Metoda siecznych. 11 1 0 + = ii i ii ii EEE i ii iiii E EE = + 1 11 )( E=0 E Ei-1 Ei i-1 i 369.6100 194.5100 11 == == ii ii E E 16.10)369.6(369.6194.5 )100100(1001 = + =+i 17.10)7.04cos(4)'7.04sin( == Przypomnienie: Wniosek: szybka zbieno ju po pierwszej iteracji. Metoda rnic skoczonych W metodzie tej, pochodne w rwnaniu rniczkowym zastpujemy rnicami skoczonymi. Dziedzin rozwiza [a,b] dzielimy na N przedziaw o dugoci h = (b-a) / N. Mamy N+1 punktw. Dla kadego z nich zapisujemy rwnanie rnicowe czyli rwnanie algebraiczne. Mamy zatem ukad rwna algebraicznych, ktry rozwizujemy jedn z metod omwionych na wykadzie nr 3. x =a1 xi xi+1xi-1 x =bN+1 h x y Przykad: 0)1(,0)0(,12 2 === yydx yd Drug pochodn mona przybliy za pomoc trjpunktowych rnic centralnych (zwykle, cho niekoniecznie). Mamy zatem Rozwizanie analityczne:22 )(2 xxxy = By otrzyma rozwizanie numeryczne najpierw dyskretyzujemy rwnanie dla punktw [x0, x1, , xN], gdzie x0 = 0, xN = 1, xi = ih. Wybierzmy h = 0.2 (N = 5). 4,3,2,1,12 2 11 ==+ + i h yyy iii dla Zwrmy uwag, e s to rwnania tylko dla punktw wewntrznych. Ukad nasz ma 4 niewiadome i 4 rwnania: . h y y- y , h y y- y , h y y-y , h y y -y 2543 2432 2321 2210 2 2 2 2 =+ =+ =+ =+ Macierz powysza jest przykadem macierzy trjdiagonalnej. Cho powyszy ukad mona rozwiza jedn z metod omwionych na wykadzie 3 (np. metod Gaussa), to szczeglna posta tej macierzy pozwala zredukowa liczb oblicze z n3 do n. = 5 2 2 2 0 2 4 3 2 1 21 121 121 12 yh h h yh y y y y = 04.0 04.0 04.0 04.0 .y y , yy y , y yy , yy 1004 43 1004 432 1004 321 1004 21 2 2 2 2 = =+ =+ =+ . y y , yy y , y y , yy 1004 43 1004 432 1006 3223 1002 221 1 2 2 = =+ =+ = . y y , yy , y y , yy 1004 43 1008 4334 1004 332 2 1002 221 1 2 = =+ = = . y , yy , yy , yy 101 445 1006 443 3 1004 332 2 1002 221 1 = = = = .)( )( )( 1008 10012 21 1002 1 10012 10012 32 1004 2 10012 1008 43 1006 3 1008 4 =+= =+= =+= = y , y , y ,y 0.0 0.2 0.4 0.6 0.8 1.0-0.14 -0.12 -0.10 -0.08 -0.06 -0.04 -0.02 0.00 y x Algorytm Thomasa dla ukadw z macierz trjdiagonaln Dany jest ukad rwna: niYycbya iiiiii ...,,2,111 ==++ + dla = nnnn n Y Y Y y y y ba c cba cba cb . . . . 0......0 ............... 0...0 0...0 0......0 2 1 2 1 1 333 222 11 lub w postaci macierzowej 0,01 == nca Algorytm skada si z dwch faz: Faza eliminacji wprzd wykonujc dla rwna od i = 1 do n eliminacj niewiadomych uzyskujemy ostatnie rwnanie (i = n) z jedn tylko niewiadom, ktr moemy wyznaczy. Faza eliminacji wstecz korzystajc z wyznaczonej w rwnaniu i+1 niewiadomej yi+1wyznaczamy z rwnania i niewiadom yi, a do otrzymania wartoci y1. ()()()() A zatem szukamy schematu postaci: iiii yy +=1 Podstawiajc ten schemat do rwnania () otrzymujemy: ( ) iiiiiiiii Yycybya =+++ +1 Wyznaczajc yi iii iiii iii ii ba aYy ba cy + + + = + 1 Porwnujc z otrzymujemy iii iiii iii ii ba aY ba c + = + = ++ 11 , Otrzymalimy rwnanie rekurencyjne na poszukiwane wspczynniki i . Wspczynniki pocztkowe 1 i 1 nie maj znaczenia, bo mnoone s przez a1=0. Musimy jeszcze zna yn, by mc rozpocz iteracyjne obliczanie niewiadomych. Podstawiajc pierwsze rwnanie schematu do rwnania () otrzymujemy: nnnn yy +=1 ( ) nnnnnnn Yybya =++ Czyli nnn nnnn ba aYy + = iii iiii ba aY + =+ 1Warto zauway, e skoro to wystarczy przyj () 01 =+ny () , by mc bezporednio skorzysta ze wzoru (). Podsumowujc: iii1-i 1n iii iii1i iii i1i 11 xy 1..2ni for 0;y ba aY, ba c 1..ni for przykad //na0; += += = + = + = = == + ++ Algorytm Thomasa jest niezawodny, gdy macierz jest diagonalnie dominujca nicab iii ...,,1=+ Przykad: radiator prtowy Rwnanie opisujce rozkad temperatury wzdu dugoci prta: TA TBTS L x LxTTkA Ph dx TdS C C = 0,0)(2 2 hc wspczynnik wnikania ciepa P obwd prta k wspczynnik przewodzenia ciepa Ac pole poprzecznego przekroju prta K293T K293T(L) K473T(0) m0.1L m101.6A W/m/K 240k m 0.016 P /KW/m40h S 25 c 2 c = = = = = = = = 0)(2 2 = STTdx Td 0)( 22 11 =+ + Siiii TT h TTT Siii ThTThT 2 1 2 1 )2( =++ + Podzielmy domen rozwiza na 5 czci (h = L / 5 = 2 cm) SThTThTi 2 32 2 1 )2(2 =++= )()2( 12 32 2TThTTh S +=++ czyli SThTThTi 2 43 2 2 )2(3 =++= SThTThTi 2 54 2 3 )2(4 =++= SThTThTi 2 65 2 4 )2(5 =++= )()2( 62 45 2TThTTh S +=++ czyli + + = + + + + )( )( )2(100 1)2(10 01)2(1 001)2( 6 2 2 2 1 2 5 4 3 2 2 2 2 2 TTh Th Th TTh T T T T h h h h S S S S Mamy zatem ukad 4 rwna z czterema niewiadomymi. 0.00 0.02 0.04 0.06 0.08 0.10 300 350 400 450 500 Rozwizujc ten ukad za pomoc algorytmu Thomasa otrzymujemy: Metoda rnic skoczonych dla nieliniowych rwna rniczkowych Najciekawsze problemy wspczesnej fizyki s nieliniowe. Nieliniowe problemy brzegowe dyskretyzujemy w podobny sposb. Wynikiem jest jednak ukad nieliniowych rwna algebraicznych. Metody rozwizywania takich rwna nieliniowych omwilimy na wykadzie nr 2. Najbardziej wydajne obliczeniowo s w tym przypadku metody iteracyjne. Istnieje jednak potencjalny problem zwizany ze zbienoci schematu iteracyjnego. Metoda punktu staego Ukad rwna nieliniowych mona zapisa w postaci ][][]][[ bya =+ [a] macierz wspczynnikw [] wektor nieliniowych wyrazw bdcych funkcj niewiadomych yi [b] - wektor znanych wielkoci staych Spord wielu sposobw konstruowania procedury iteracyjnej wybierzmy najprostszy ][][]][[ bya =+ kkbya ][][]][[ 1 =+ obliczone na podstawie wczeniejszego kroku k Jeli liczba punktw jest maa moemy macierz [a] odwrci. Jeli nie, moemy skorzysta z metody eliminacji Gaussa albo Thomasa (dla macierzy trjdiagonalnej) ( )kk bay ][][][][ 11 = + Przykad: radiator prtowy Gdy uwzgldnimy wyraz odpowiedzialny na wypromieniowywanie ciepa, rwnanie opisujce rozkad temperatury wzdu dugoci prta uzyska posta: LxTTkA PTT kA Ph dx TdS C S C C = 0,0)()( 442 2 - wzgldna zdolno emisyjna - staa Stefana-Boltzmanna 0)()(2 44 2 11 =+ + SiBSiAiii TTTT h TTT Rwnanie zdyskretyzowane: )()2( 421422 1 SBSAiiBiAi TThTThThT +=++ + = + + + + + )( )( )( )( )2(100 1)2(10 01)2(1 001)2( 4 5 2 4 4 2 4 3 2 4 2 2 5 4 3 2 2 2 2 2 Th Th Th Th T T T T h h h h B B B B A A A A + + + + = 6 42 42 42 1 42 )( )( )( )( TTTh TTh TTh TTTh SBSA SBSA SBSA SBSA ( )kk baT ][][][][ 11 = +Teraz stosujemy procedur iteracyjn: 0.00 0.02 0.04 0.06 0.08 0.10 280 320 360 400 440 480 T [K ] x [m] krok 0 krok 1 krok 2 krok 3 0 0.02 0.04 0.06 0.08 0.1 krok 0 473 400 400 400 400 293 krok 1 473 423.2293 382.8297 349.1078 319.8155 293 krok 2 473 423.3492 383.3225 349.8507 320.4519 293 krok 3 473 423.344 383.3132 349.8409 320.4456 293 Zakadamy hx = hy = h [siatka kwadratowa] u(xi,yj) = ui,j u(xi+h, yj+h) = ui+1,j+1, xi xi+1xi-1 yjyj-1 yj+1 ( ) [ ]jijijixji uuuhyxu ,1,,12 21 ,'' + + 02 2 2 2 = + y u x uui,j hx hy Dwuwymiarowe zagadnienie brzeg

Recommended