Baixe o app para aproveitar ainda mais
Prévia do material em texto
O Me´todo das Diferenc¸as Finitas Equac¸o˜es na˜o lineares (Navier-Stokes) Prof: Rafael Gabler Gontijo November 29, 2014 Abstract Esse documento apresenta um esquema utilizado para soluc¸a˜o de sis- temas de equac¸o˜es na˜o lineares. O objetivo e´ que o aluno possa uti- lizar o Me´todo de Newton para resolver um sistema de equac¸o˜es comum na Mecaˆnica dos Fluidos, associado a` forma discretizada da equac¸a˜o de Navier-Stokes em um problema de camada limite. O esquema se baseara´ em um modelo dida´tico de uma malha mais simples contendo 16 no´s a t´ıtulo de ilustrac¸a˜o. 1 Introduc¸a˜o e motivac¸a˜o Em muitos problemas da f´ısica precisamos resolver equac¸o˜es diferenciais or- dina´rias ou parciais para entendermos como determinada grandeza se com- porta. Muitas leis fundamentais da f´ısica cla´ssica sa˜o expressas em linguagem matema´tica atrave´s de equac¸o˜es diferenciais. Essas equac¸o˜es geralmente esta˜o associadas ao comportamento de func¸o˜es escalares, vetoriais ou tensoriais, no espac¸o e no tempo. Nesse sentido o conceito de campo e´ muito utilizado para descrever o comportamento de um sistema f´ısico. Quando um determinado sistema f´ısico e´ descrito em linguagem matema´tica atrave´s de equac¸o˜es diferenciais que possuem soluc¸a˜o anal´ıtica, esta e´ a te´cnica utilizada para entender o comportamento deste. A vantagem de uma soluc¸a˜o anal´ıtica de uma equac¸a˜o diferencial e´ que esta e´ exata e fornece uma expressa˜o matema´tica como resposta. Esta expressa˜o por sua vez pode ser interpretada para gerar uma teoria f´ısica. E´ um consenso entre os principais cientistas de que a soluc¸a˜o anal´ıtica e´ sem du´vidas a busca primeira e primordial para o desenvolvimento de uma teoria cient´ıfica. Entrentanto, frequentemente nos de- paramos com problemas regidos por equac¸o˜es diferenciais com elevado grau de complexidade que impedem o desenvolvimento de uma soluc¸a˜o anal´ıtica. Nesses casos podemos abordar o problema com dois focos poss´ıveis: uma soluc¸a˜o nume´rica (computacional) ou um estudo experimental para buscar correlac¸o˜es emp´ıricas que representem o comportamento do sistema em estudo. A vantagem de uma abordagem experimental consiste no fato de que se o experimento for conduzido de forma meticulosa e com o menor grau poss´ıvel de interfereˆncia do sistema de medic¸a˜o no sistema, ela representa uma exce- lente visa˜o do que poder´ıamos chamar de realidade. As desvantagens de uma abordagem experimental sa˜o: necessidade de ana´lises estat´ısticas custosas para 1 garantir a repetibilidade dos resultados, alto custo (dependendo da grandeza f´ısica a ser mensurada), interfereˆncia do dispositivo de medida (senores, sondas) no sistema, em muitos casos as medic¸o˜es sa˜o indiretas e dependem da inter- pretac¸a˜o do experimentalista, dificuldade de controle de fatores externos que podem influenciar no experimento (som, luz, temperatura, pressa˜o, ru´ıdos da rede ele´trica, impurezas, entre outros). A abordagem nume´rica muitas vezes e´ o caminho preferido por muitos ci- entistas (ainda mais com a enorme evoluc¸a˜o tecnlo´gica recente que impulsiona o poder de ca´lculo de nossos computadores modernos) por apresentar vanta- gens como: facilidade no controle dos paraˆmetros f´ısicos do problema, acesso facilitado a inu´meras ferramentas de pre´ e po´s processamento de dados, di- versidade de linguagens de programac¸a˜o, acesso amplo a` inu´meras te´cnicas de discretizac¸a˜o espacial para soluc¸a˜o de problemas de valor de contorno, entre outros. Evidamentemente o me´todo nume´rico tambe´m apresenta desvantagens, como: depedeˆncia dos resultados de acordo com a resoluc¸a˜o espacial e temporal adotada, necessidade eventual de calibrac¸a˜o de paraˆmetros em aproximac¸o˜es de func¸o˜es cont´ınuas por valores discretos e (provavelmente a principal delas) a impossibilidade de obtenc¸a˜o de uma expressa˜o matema´tica geral como soluc¸a˜o da equac¸a˜o diferencial a ser resolvida. Tudo que o me´todo nume´rico (como o pro´prio nome diz) fornece sa˜o nu´meros. Nu´meros estes que podem ser or- ganizados na forma de tabelas, gra´ficos e imagens de campo para estimar o comportamento f´ısico aproximado de determinado sistema. Figure 1: Exemplo do uso de um me´todo computacional para a soluc¸a˜o da equac¸a˜o de Navier-Stokes para a determinac¸a˜o do campo de velocidades de um escoamento na sa´ıda de uma turbina. De um modo geral a lo´gica principal por tra´s das diversas estrate´gias de soluc¸a˜o de equac¸o˜es diferenciais cont´ınuas atrave´s de um me´todo nume´rico con- siste em discretizar um domı´nio espacial de ca´lculo cont´ınuo (com infinitos pon- tos) em um domı´nio espacial discreto (com um nu´mero finito de pontos). Essa transformac¸a˜o tambe´m se aplica ao tempo (em problemas transientes) que em 2 sistemas reais flui de forma cont´ınua, mas em sistemas nume´ricos deve ser dis- cretizado em diferentes instantes sucessivos descont´ınuos. Essa transformac¸a˜o cont´ınuo-discreto e´ tambe´m responsa´vel pela trans- formac¸a˜o da necessidade de soluc¸a˜o de uma u´nica equac¸a˜o diferencial em N equac¸o˜es alge´bricas para cada instante de tempo (aqui N representa o nu´mero de no´s da malha que representa o espac¸o da soluc¸a˜o discretizado). Cada equac¸a˜o alge´brica se refere a` maneira com a qual o valor de uma determinada varia´vel em um no´ i qualquer da malha de ca´lculo se relaciona com seus no´s vizinhos. Essa relac¸a˜o evidentemente e´ dada pela discretizac¸a˜o da equac¸a˜o diferencial gover- nante do problema, que por sua vez e´ deduzida com base em um (ou um conjunto de) princ´ıpio(s) f´ısico(s). Como um no´ de refereˆncia qualquer da malha possui ao menos 4 vizinhos (em sistemas bidimensionais), este no´ em uma ana´lise pos- terior sera´ o vizinho de um outro no´ de refereˆncia. Por este motivo a conexa˜o entre os no´s da malha de ca´lculo e´ intricada, de modo que um no´ depende do outro e vice-versa. Por esse motivo esse conjunto de N equac¸o˜es discretas (resul- tado da transformac¸a˜o cont´ınuo-discreto da equac¸a˜o diferencial governante) e´ comumente expresso atrave´s de um sistema linear de equac¸o˜es. Um exemplo cla´ssico e´ o da equac¸a˜o da conduc¸a˜o de calor bidimensional permanente. Para esta equac¸a˜o temos a seguinte expressa˜o cont´ınua: ∇2T = 0 → ∂2T ∂x2 + ∂2T ∂y2 = 0, eq. cont´ınua, (1) em que ∇ e´ o operador nabla, T e´ o campo de temperatura, x e y representam coordenadas espaciais Eulerianas. Para uma malha na qual o espac¸amento entre os no´s na direc¸a˜o x, dado por ∆x e´ igual ao espac¸amento y dado por ∆y, a versa˜o discretizada (via me´todo das diferenc¸as finitas) da equac¸a˜o (1) e´ dada por Ti − 1 4 (Ti−1 + Ti+1 + Ti−n + Ti+n) = 0 (2) aqui n se refere ao nu´mero de no´s na direc¸a˜o x. A equac¸a˜o (2) pode ser escrita na forma de um sistema linear para os no´s internos como · · · ... ... ... ... ... ... ... · · · · · · −0.25 · · · −0.25 1 −0.25 · · · −0.25 · · · · · · ... ... ... ... ... ... ... · · · · ... Ti−n ... Ti−1 Ti Ti+1 ... Ti+n ... = ... ... ... ... 0 ... ... ... ... (3) O preenchimento das outras linhas da matriz A (primeiro termo do lado es- querdo), do vetor x (segundo termo do lado esquerdo) e do vetor B (termo do lado direito da equac¸a˜o) dependem das condic¸o˜es de contorno do problema. 3 De qualquer forma a soluc¸a˜o da equac¸a˜o cont´ınua (1) em um domı´nio dis- cretizado com N no´s passa pela soluc¸a˜o de um sistema linear de ordem N do tipo A·x = B. Existem diversos me´todos dispon´ıveis para a soluc¸a˜o de sistemas lineares. Alguns desses me´todos demandam a inversa˜o da matriz A (proced- imento muito caro do ponto de vista computacional), outrosexigem soluc¸o˜es iterativas. O problema da soluc¸a˜o de um sistema linear por me´todos convencionais para soluc¸a˜o destes tipos de sistema e´ que estas sa˜o va´lidas para a soluc¸a˜o de sistemas de equac¸o˜es discretizadas oriundos de equac¸o˜es diferenciais cont´ınuas lineares. Em outras palavras: na˜o sa˜o va´lidas para sistemas na˜o lineares. Algu- mas equac¸o˜es diferenciais importantes que regem uma gama enorme de sistemas f´ısicos sa˜o na˜o lineares. Um exemplo cla´ssico e´ a equac¸a˜o de Navier-Stokes que, para um fluido incompress´ıvel, e´ dada por ρ ( ∂u ∂t + u · ∇u ) = −∇p+ η∇2u+ ρg, (4) em que ρ e η representam a massa espec´ıfica e a viscosidade dinaˆmica do fluido respectivamente, p e u denotam os campos Eulerianos de pressa˜o e velocidade do escoamento, t e´ a varia´vel tempo e g representa o vetor acelerac¸a˜o do campo gravitacional. Note que o termo u ·∇u torna a equac¸a˜o (4) uma EDP na˜o linear de segunda ordem na˜o homogeˆnea. O elevado grau de complexidade da equac¸a˜o (4) nos impossibilita uma soluc¸a˜o anal´ıtica geral, de modo que esta so´ apresenta soluc¸a˜o em casos assinto´ticos muito particulares. Dessa forma muitas vezes a Mecaˆnica dos Fluidos avanc¸a em termos de compreensa˜o f´ısica atrave´s de sim- ulac¸o˜es computacionais baseadas na soluc¸a˜o da equac¸a˜o (4) com condic¸o˜es ini- ciais e de contorno bem definidas (problema bem posto) via me´todos nume´ricos de discretizac¸a˜o temporal. O grande inconveniente de discretizarmos a equac¸a˜o (4) cheia (com todos os termos) e´ o surgimento de um sistema de equac¸o˜es na˜o lineares. Este documento apresenta um esquema dida´tico mostrando como pode- mos resolver uma equac¸a˜o diferencial parcial na˜o linear atrave´s do me´todo das diferenc¸as finitas combinado com te´cnicas de soluc¸a˜o de sistemas de equac¸o˜es na˜o lineares. 2 O me´todo de Newton - ra´ızes de func¸o˜es na˜o lineares Existem diversos me´todos utilizados para a soluc¸a˜o de equac¸o˜es na˜o lineares, que va˜o desde uso de algoritmos gene´ticos, me´todo da descida ı´ngreme, me´todos de ponto fixo, entre outros. Um me´todo de fa´cil compreensa˜o e ao mesmo tempo muito eficaz consiste no me´todo de Newton. Geralmente alunos de cursos de ca´lculo nume´rico esta˜o familiarizados com o me´todo de Newton para a deter- minac¸a˜o de zeros de func¸o˜es e em certo sentido e´ essa a ideia por tra´s de seu uso para resolver sistemas de equac¸o˜es na˜o lineares. A ideia do me´todo de Newton e´ que uma func¸a˜o f(x) qualquer (linear ou na˜o) pode ser aproximada por uma se´rie de Taylor no ponto x = x0 +∆x por seu valor e de sua derivada no ponto x0 como f(x) ≈ f(x0) + df dx (x0)∆x (5) 4 aqui ∆x = x − x0. Caso desejemos encontrar alguma ra´ız da func¸a˜o f(x), ou seja, determinar o valor (ou valores) de x que fazem com que f(x) = 0, igualamos a equac¸a˜o (5) e isolamos x, este procedimento fornece x = x0 − f(x0) f ′(x0) (6) aqui f ′ (x0) denota a primeira derivada de f com respeito a x avaliada em x0. Como o me´todo se baseia na aproximac¸a˜o dada por (5) o uso da equac¸a˜o (6) na˜o fornecera´ a resposta exata com apenas uma u´nica computac¸a˜o. Na realidade o me´todo de Newton e´ um me´todo iterativo, em que xi+1 = xi − f(xi) f ′(xi) , i = 0, 1, 2, · · · (7) A partir do momento em que a diferenc¸a entre dois valores consecutivos de x for suficientemente pequena (dentro de uma toleraˆncia pre´-estabelecida) dizemos que o me´todo convergiu. Considere por exemplo que desejemos determinar uma ra´ız da func¸a˜o dada por f(x) = x2 − 4, (8) aplicando o me´todo descrito em (7) para a equac¸a˜o (8) obtemos xi+1 = xi − xi 2 − 4 2xi , i = 0, 1, 2, · · · (9) Sabemos que a equac¸a˜o (8) possui duas ra´ızes reais, dadas por 2 e -2. Con- siderando x0 = 0.5 podemos montar a seguinte tabela: i xi xi+1 |xi+1 − xi| 0 0.5 4.25 3.75 1 4.25 2.5956 1.6544 2 2.5956 2.0683 0.5273 3 2.0683 2.0011 0.0672 4 2.0011 2.0000003183 0.0011284441 5 2.0000003183 2.0 3.18 ×10−7 Table 1: Determinac¸a˜o de uma das ra´ızes de uma equac¸a˜o na˜o linear alge´brica de segundo grau simples via Me´todo de Newton. Note que a partir da sexta iterac¸a˜o a diferenc¸a entre dois valores consecutivos da ra´ız obtida e´ da ordem de 10−7. Essa tabela mostra a efica´cia do me´todo e sua ra´pida convergeˆncia. O uso do me´todo de Newton para a soluc¸a˜o de sistemas na˜o lineares se baseia na generalizac¸a˜o deste para encontrar N ra´ızes de N equac¸o˜es na˜o lineares organizadas na forma de um sistema na˜o linear. Um sistema na˜o linear de N equac¸a˜o acopladas pode ser organizado em termos vetoriais como uma func¸a˜o f(x) = 0. Essa func¸a˜o f e´ um vetor com N componentes em que N e´ o nu´mero de equac¸o˜es a ser resolvida. Como exemplo considere o seguinte sistema de duas equac¸o˜es na˜o lineares f(x) = ( f1 f2 ) = ( x1 + 2x2 − 2 x21 + 4x 2 2 ) = ( 0 0 ) , (10) 5 considere ainda a confecc¸a˜o da seguinte matriz J J = ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2 = 1 2 2x1 8x2 . (11) Essa matriz e´ denominada matriz Jacobiana do sistema e e´ utilizada para escr- ever a equac¸a˜o de recorreˆncia (9) em forma vetorial como xi+1 = xi − J(xi) −1 · f(xi). (12) Note que na equac¸a˜o (12) surge a matriz inversa de J . Isso na˜o significa necessariamente que precisamos inverter a matriz J para resolver o sistema. Chamando o termo −J(xi) −1 · f(xi) de si temos si = −J(xi) −1 · f(xi) (13) multiplicando os dois lados da equac¸a˜o (13) por J(xi) obtemos J(xi) · si = −J(xi) · J(xi) −1 · f(xi) = −I · f(xi) = −f(xi), (14) dessa forma temos xi+1 = xi + si, i = 0, 1, 2, · · · (15) com si dado pela soluc¸a˜o do sistema linear J(xi) · si = −f(xi). (16) Para o exemplo em questa˜o considere um vetor inicial x0 = [1 2]T , dessa forma f(x0) = [3 13]T e J(x0) = 1 2 2 16 (17) Resolvendo o sistema 1 2 2 16 · s0 = ( −3 13 ) , (18) obtemos s0 = [−1.83 −0.58]T e x1 = [−0.83 1.42]T . Para um pro´ximo ponto temos f(x1) = [0 4.72]T e J(x1) = 1 2 −1.67 11.3 (19) Resolvendo o sistema 1 2 −1.67 11.3 · s1 = ( 0 4.72 ) , (20) obtemos s1 = [0.64 − 0.32]T e x1 = [−0.19 1.1]T , repetindo o procedimento diversas vezes a soluc¸a˜o convirgira´ para x = [0 1]T . 6 Note que o uso do Me´todo de Newton para a soluc¸a˜o de sistemas na˜o lineares consiste na transformac¸a˜o de um u´nico sistema na˜o linear em diversos sitemas lineares resolvidos de modo progressivo em um procedimento iterativo ate´ a convergeˆncia. A ideia por tra´s deste me´todo e´ a de que uma func¸a˜o na˜o linear pode ser aproximada como uma soma de diferentes func¸o˜es lineares em torno de pequenos intervalos. A pro´xima sec¸a˜o mostra um exemplo do uso deste me´todo na soluc¸a˜o da equac¸a˜o de Navier-Stokes para um problema de camada limite laminar sobre placa plana de um fluido incompress´ıvel. 3 Exemplo - Soluc¸a˜o de uma EDP na˜o linear Considere a seguinte equac¸a˜o diferencial parcial u ∂u ∂x + v ∂u ∂x = 1 Re ∂2u ∂y2 (21) a equac¸a˜o (21) e´ conhecida como equac¸a˜o de Navier-Stokes (sua componente apenas na direc¸a˜o x para um problema de camada limite). Considere a seguinte imagem ilustrativa que representa o domı´nio de ca´lculo discretizado no qual desejamos resolver a equac¸a˜o (21). A figura (3) ilustra um modelo esquema´tico ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ���������������������������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� ����������������� 1 2 3 4 5 6 7 8 9 10 11 12 14 15 1613 u = 1.0 u = 0.0 u = 1.0 x y i i+1i−1 i+n i−n du/dx = 0 Figure 2: O domı´nio de ca´lculo discretizado consiste em uma malha com 16 no´s, sendo 4 internos e 12 externos sujeitos a`s condic¸o˜es de contorno descritas na figura. do problema de camda limite a ser abordado. As func¸o˜es escalares u e v sa˜o func¸o˜es de x e y (coordenadas espaciais), de modo que u = u(x, y) e v = v(x, y) (22) A versa˜o discretizada da equac¸a˜o (21) e´ dada por ui ( ui − ui−1 ) ∆x + vi ( ui − ui−n ) ∆x = 1 Re ( ui+n + ui−n − 2ui ) ∆y2 (23) 7 Figure 3: Modelo esquema´tico do problema de camada limite hidrodinaˆmica. isolando ui temos ui ( ui ∆x + vi ∆y + 2 Re∆y2 − ui−1 ∆x ) − viui−n ∆y − ui+n Re∆y2 − ui−n Re∆y2 = 0 (24) Note que a equac¸a˜o acima se refere a` um no´ i qualquer e a maneira com a qual este se relaciona com seus vizinhos. Para fins de simplificac¸a˜o considere o caso em que v = 0, de modo que ui ( ui ∆x + 2 Re∆y2 − ui−1 ∆x ) − ui+n Re∆y2 − ui−n Re∆y2 = 0 (25) Podemos escrever esse sistema de equac¸o˜es na˜o lineares para a malha em questa˜o (com 16 no´s) da forma 8 f(u) = u1 u2 u3 u4 u5 − 1 u2 6 ∆x − u6u5 ∆x + 2u6 Re∆y2 − u10 Re∆y2 − u2 Re∆y2 u2 7 ∆x − u7u6 ∆x + 2u7 Re∆y2 − u11 Re∆y2 − u3 Re∆y2 u8 − u7 u9 − 1 u2 10 ∆x − u10u9 ∆x + 2u10 Re∆y2 − u14 Re∆y2 − u6 Re∆y2 u2 11 ∆x − u11u10 ∆x + 2u11 Re∆y2 − u15 Re∆y2 − u7 Re∆y2 u12 − u11 u13 − 1 u14 − 1 u15 − 1 u16 − 1 = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (26) Considere a matriz Jacobiana J do sistema definida como 9 J = ∂f1 ∂u1 ∂f1 ∂u2 ∂f1 ∂u3 ∂f1 ∂u4 ∂f1 ∂u5 ∂f1 ∂u6 ∂f1 ∂u7 ∂f1 ∂u8 ∂f1 ∂u9 ∂f1 ∂u10 ∂f1 ∂u11 ∂f1 ∂u12 ∂f1 ∂u13 ∂f1 ∂u14 ∂f1 ∂u15 ∂f1 ∂u16 ∂f2 ∂u1 ∂f2 ∂u2 ∂f2 ∂u3 ∂f2 ∂u4 ∂f2 ∂u5 ∂f2 ∂u6 ∂f2 ∂u7 ∂f2 ∂u8 ∂f2 ∂u9 ∂f2 ∂u10 ∂f2 ∂u11 ∂f2 ∂u12 ∂f2 ∂u13 ∂f2 ∂u14 ∂f2 ∂u15 ∂f2 ∂u16 ∂f3 ∂u1 ∂f3 ∂u2 ∂f3 ∂u3 ∂f3 ∂u4 ∂f3 ∂u5 ∂f3 ∂u6 ∂f3 ∂u7 ∂f3 ∂u8 ∂f3 ∂u9 ∂f3 ∂u10 ∂f3 ∂u11 ∂f3 ∂u12 ∂f3 ∂u13 ∂f3 ∂u14 ∂f3 ∂u15 ∂f3 ∂u16 ∂f4 ∂u1 ∂f4 ∂u2 ∂f4 ∂u3 ∂f4 ∂u4 ∂f4 ∂u5 ∂f4 ∂u6 ∂f4 ∂u7 ∂f4 ∂u8 ∂f4 ∂u9 ∂f4 ∂u10 ∂f4 ∂u11 ∂f4 ∂u12 ∂f4 ∂u13 ∂f4 ∂u14 ∂f4 ∂u15 ∂f4 ∂u16 ∂f5 ∂u1 ∂f5 ∂u2 ∂f5 ∂u3 ∂f5 ∂u4 ∂f5 ∂u5 ∂f5 ∂u6 ∂f5 ∂u7 ∂f5 ∂u8 ∂f5 ∂u9 ∂f5 ∂u10 ∂f5 ∂u11 ∂f5 ∂u12 ∂f5 ∂u13 ∂f5 ∂u14 ∂f5 ∂u15 ∂f5 ∂u16 ∂f6 ∂u1 ∂f6 ∂u2 ∂f6 ∂u3 ∂f6 ∂u4 ∂f6 ∂u5 ∂f6 ∂u6 ∂f6 ∂u7 ∂f6 ∂u8 ∂f6 ∂u9 ∂f6 ∂u10 ∂f6 ∂u11 ∂f6 ∂u12 ∂f6 ∂u13 ∂f6 ∂u14 ∂f6 ∂u15 ∂f6 ∂u16 ∂f7 ∂u1 ∂f7 ∂u2 ∂f7 ∂u3 ∂f7 ∂u4 ∂f7 ∂u5 ∂f7 ∂u6 ∂f7 ∂u7 ∂f7 ∂u8 ∂f7 ∂u9 ∂f7 ∂u10 ∂f7 ∂u11 ∂f7 ∂u12 ∂f7 ∂u13 ∂f7 ∂u14 ∂f7 ∂u15 ∂f7 ∂u16 ∂f8 ∂u1 ∂f8 ∂u2 ∂f8 ∂u3 ∂f8 ∂u4 ∂f8 ∂u5 ∂f8 ∂u6 ∂f8 ∂u7 ∂f8 ∂u8 ∂f8 ∂u9 ∂f8 ∂u10 ∂f8 ∂u11 ∂f8 ∂u12 ∂f8 ∂u13 ∂f8 ∂u14 ∂f8 ∂u15 ∂f8 ∂u16 ∂f9 ∂u1 ∂f9 ∂u2 ∂f9 ∂u3 ∂f9 ∂u4 ∂f9 ∂u5 ∂f9 ∂u6 ∂f9 ∂u7 ∂f9 ∂u8 ∂f9 ∂u9 ∂f9 ∂u10 ∂f9 ∂u11 ∂f9 ∂u12 ∂f9 ∂u13 ∂f9 ∂u14 ∂f9 ∂u15 ∂f9 ∂u16 ∂f10 ∂u1 ∂f10 ∂u2 ∂f10 ∂u3 ∂f10 ∂u4 ∂f10 ∂u5 ∂f10 ∂u6 ∂f10 ∂u7 ∂f10 ∂u8 ∂f10 ∂u9 ∂f10 ∂u10 ∂f10 ∂u11 ∂f10 ∂u12 ∂f10 ∂u13 ∂f10 ∂u14 ∂f10 ∂u15 ∂f10 ∂u16 ∂f11 ∂u1 ∂f11 ∂u2 ∂f11 ∂u3 ∂f11 ∂u4 ∂f11 ∂u5 ∂f11 ∂u6 ∂f11 ∂u7 ∂f11 ∂u8 ∂f11 ∂u9 ∂f11 ∂u10 ∂f11 ∂u11 ∂f11 ∂u12 ∂f11 ∂u13 ∂f11 ∂u14 ∂f11 ∂u15 ∂f11 ∂u16 ∂f12 ∂u1 ∂f12 ∂u2 ∂f12 ∂u3 ∂f12 ∂u4 ∂f12 ∂u5 ∂f12 ∂u6 ∂f12 ∂u7 ∂f12 ∂u8 ∂f12 ∂u9 ∂f12 ∂u10 ∂f12 ∂u11 ∂f12 ∂u12 ∂f12 ∂u13 ∂f12 ∂u14 ∂f12 ∂u15 ∂f12 ∂u16 ∂f13 ∂u1 ∂f13 ∂u2 ∂f13 ∂u3 ∂f13 ∂u4 ∂f13 ∂u5 ∂f13 ∂u6 ∂f13 ∂u7 ∂f13 ∂u8 ∂f13 ∂u9 ∂f13 ∂u10 ∂f13 ∂u11 ∂f13 ∂u12 ∂f13 ∂u13 ∂f13 ∂u14 ∂f13 ∂u15 ∂f13 ∂u16 ∂f14 ∂u1 ∂f14 ∂u2 ∂f14 ∂u3 ∂f14 ∂u4 ∂f14 ∂u5 ∂f14 ∂u6 ∂f14 ∂u7 ∂f14 ∂u8 ∂f14 ∂u9 ∂f14 ∂u10 ∂f14 ∂u11 ∂f14 ∂u12 ∂f14 ∂u13 ∂f14 ∂u14 ∂f14 ∂u15 ∂f14 ∂u16 ∂f15 ∂u1 ∂f15 ∂u2 ∂f15 ∂u3 ∂f15 ∂u4 ∂f15 ∂u5 ∂f15 ∂u6 ∂f15 ∂u7 ∂f15 ∂u8 ∂f15 ∂u9 ∂f15 ∂u10 ∂f15 ∂u11 ∂f15 ∂u12 ∂f15 ∂u13 ∂f15 ∂u14 ∂f15 ∂u15 ∂f15 ∂u16 ∂f16 ∂u1 ∂f16 ∂u2 ∂f16 ∂u3 ∂f16 ∂u4 ∂f16 ∂u5 ∂f16 ∂u6 ∂f16 ∂u7 ∂f16 ∂u8 ∂f16 ∂u9 ∂f16 ∂u10 ∂f16 ∂u11 ∂f16 ∂u12 ∂f16 ∂u13 ∂f16 ∂u14 ∂f16 ∂u15 ∂f16 ∂u16 (27) Determinando as derivadas especificadas em (27) para o sistema definido em (26) obtemos 10 J = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 [ − 1 Re∆y2 ] 0 0 [ − u6 ∆x ] [ 1 ∆x (2u6 − u5) + 2 Re∆y2 ] 0 0 0 [ − 1 Re∆y2 ] 0 0 0 0 0 0 0 0 [ − 1 Re∆y2 ] 0 0 [ − u7 ∆x ] [ 1 ∆x (2u7 − u6) + 2 Re∆y2 ] 0 0 0 [ − 1 Re∆y2 ] 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 (28) 1 1 Considere o caso em que ∆x = ∆y = 0.1 e Re = 100, essa matriz e´ dada por J = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −10u6 10 (2u6 − u5) + 2 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −10u7 10 (2u7 − u6) + 2 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 (29) O me´todo de soluc¸a˜o consiste em supor uma condic¸a˜o inicial u0 para o vetor u, substituir essa condic¸a˜o na matriz J(u0) e resolver o seguinte sistema linear J(u0) · u ∗ 0 = −f(u0) (30) em seguida no´s atualizamos o novo valor do vetor u1 fazendo u1 = u0 + u ∗ 0 , (31) com esse novo valor de u1 no´s resolvemos o seguinte sistema J(u1) · u ∗ 1 = −f(u1), (32) em seguida fazemos u2 = u1 + u ∗ 1 , (33) e assim sucessivamente ate´ que o me´todo convirja para um valor dentro de umadeterminada faixa de toleraˆncia. Podemos resumir o me´todo nos seguintes passos 1. Supor um valor inicial para o vetor u0 2. Determinar a matriz Jacobiana J(u0) 3. Determinar o vetor f (u0) 4. Resolver o sistema linear J(u0) · u ∗ 0 = −f(u0) 12 5. Corrigir o pro´ximo valor de u1 com u1 = u0 + u ∗ 0 6. Fazer para j de 1 a` N: • Resolver o sistema linear J(uj) · u ∗ j = −f(uj) • Corrigir o pro´ximo valor de uj com uj = uj + u ∗ j 7. Repetir o procedimento ate´ convergir; Para fins de exemplificac¸a˜o considere o seguinte valor inicial para o vetor u0 u0 = 0 0 0 0 1 1 1 1 1 1 1 1 1 1 1 1 (34) Para esses valores iniciais, com ∆x = ∆y = 0.1 e Re = 100, temos que J(u0) e´ dada por J = 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 (35) 13 e o vetor f (u0) e´ dado por f (u0) = 0 0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 (36) Dessa forma precisamos resolver inicialmente o seguinte sistema 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 −1 0 0 −10 12 0 0 0 −1 0 0 0 0 0 0 0 0 0 0 0 −1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 · u∗0 1 u∗0 2 u∗0 3 u∗0 4 u∗0 5 u∗0 6 u∗0 7 u∗0 8 u∗0 9 u∗0 10 u∗0 11 u∗0 12 u∗0 13 u∗0 14 u∗0 15 u∗0 16 = 0 0 0 0 0 −1 −1 0 0 0 0 0 0 0 0 0 (37) A soluc¸a˜o desse sistema linear pelo me´todo de Jacobi (me´todo iterativo) com 14 500 iterac¸o˜es fornece u∗ 0 = 0 0 0 0 0 −0.083916068 −0.1548242 −0.1548242 0 −0.006993004 −0.018729515 −0.018729515 0 0 0 0 (38) Fac¸amos agora a seguinte correc¸a˜o u1 = u0 + u ∗ 0 , (39) de modo que u1 = 0 0 0 0 1 0.9160839 0.8451759 0.8451759 1 0.9930070 0.9812705 0.9812705 1 1 1 1 (40) Agora precisamos resolver o seguinte sistema J(u1) · u ∗ 1 = −f(u1), (41) note que a matriz J muda a cada nova iterac¸a˜o. Montando o sistema e resol- vendo para u∗ 1 pelo mesmo me´todo (Jacobi com 500 iterac¸o˜es) obtemos o novo 15 valor de u2 dado por u2 = 0 0 0 0 1 0.9092013 0.8277108 0.8277108 1 0.9923854 0.9790677 0.9790677 1 1 1 1 (42) Repetindo esse procedimento 10 vezes obtemos u10, comparando o valor de u10 com o valor de u9 temos u9 = 0 0 0 0 1 0.9091543 0.8274714 0.8274714 1 0.9923811 0.9790405 0.9790405 1 1 1 1 u10 = 0 0 0 0 1 0.9091544 0.8274716 0.8274716 1 0.9923812 0.9790406 0.9790406 1 1 1 1 (43) 16 Para ana´lise da convergeˆncia do processo repare no vetor dado por |u10 − u9| |u10 − u9| = 0 0 0 0 0 1.1920929E− 07 1.1920929E− 07 1.1920929E− 07 0 5.9604645E− 08 5.9604645E− 08 5.9604645E− 08 0 0 0 0 (44) O valor ma´ximo dessa diferenc¸a, nesse caso 1.1920929E− 07 pode ser utilizado como o paraˆmetro que devera´ ser menor que determinada toleraˆncia pre´-definida pelo usua´rio para garantir a convergeˆncia do me´todo. 4 Bibliografia • Notas de aula do curso de fenoˆmenos de transporte - Prof. Rafael Gabler Gontijo • Micromecaˆnica e Microhidrodinaˆmica de Suspenso˜es Magne´ticas - Tese de Doutorado - 2013 - Prof. Rafael Gabler Gontijo 17 Introdução e motivação O método de Newton - raízes de funções não lineares Exemplo - Solução de uma EDP não linear Bibliografia
Compartilhar