Modelagem-Simulação-Processos-Diferencas-Finitas

Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda Cada um dos termos das derivadas parciais de segunda ordem são então substituídos por suas aproximações em diferenças centrais para obter 1Δx²ui1jk 2uijk ui1jk 1Δy²uij1k 2uijk uij1k 44 A qual pode ser rearranjada para 21Δx² 1Δy²uij 1Δx²ui1j 1Δx²ui1j 1Δy²uij1 1Δy²uij1 0 45 Essa é uma equação algébrica linear envolvendo o valor da variável dependente em cinco pontos adjacentes da malha Um objeto retangular dividido em p segmentos na direção x e q segmentos na direção y tem p 1 q 1 pontos na malha e p 1 q 1 pontos internos A Equação 45 escrita para cada uma dos pontos internos constitui um conjunto de p 1 q 1 equações lineares algébricas simultâneas com p1 q 1 4 incógnitas os quatro pontos dos cantos não aparecem nessas equações As condições de contorno fornecem informações adicionais para a solução do problema Se as condições de contorno são do tipo Dirichlet os valores da variável dependente são conhecidos em todos os pontos externos da malha Por outro lado se as condições de contorno são do tipo Neumann ou Robbins as quais especificam as derivadas parciais nas fronteiras essas condições devem ser também substituídas por aproximações por diferenças finitas Podese demonstrar isso especificandose uma condição de Neumann na fronteira esquerda isto é ux β em x 0 e todo y 46 onde β é uma constante Substituindo a derivada parcial na Equação 46 com uma aproximação da diferença central obtémse 12Δx ui1j ui1j β 47 Isso só é valido em x 0 onde i 0 portanto a Equação 47 se torna u1j u1j 2βΔx 48 Os pontos 1j estão localizados fora do objeto em questão dessa maneira u1j possuem valores fictícios Seu cálculo no entanto é necessário para a avaliação da condição de contorno de Neumann A Equação 48 escrita para todo y j 01q fornece q 1 equações adicionais mas ao mesmo tempo introduz q 1 variáveis adicionais Para balancear essa situação a Equação 45 é também escrita para q 1 pontos ao longo dessa fronteira em x 0 fornecendo o número necessário de equações independentes para a solução do problema Substituindo a derivada parcial na Equação 46 com uma diferença progressiva não requer o uso de pontos fictícios No entanto é importante utilizar a fórmula da diferença progressiva com a mesma acurácia das outras equações 54 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda 12Δxu2j 41j 3u0j β 49 ou 3u0j 4u1j u2j 2βΔx 410 A Equação 410 fornece q 1 equações adicionais sem introduzir variáveis adicionais No caso de uma condição de Robbins na fronteira esquerda na forma ux β γu em x 0 e todo y 411 onde β e γ são constantes uma derivação similar como a anterior mostra que a seguinte equação deve ser utilizada na fronteira 3 2γΔxu0j 4u1j u2j 2βΔx 412 A Equação 45 e a condição de contorno apropriada constituem um conjunto de equações algébricas lineares então os métodos de Gauss para a solução dessas equações podem ser usados A Equação 45 é um sistema predominantemente diagonal dessa forma o método de GaussSiedel é especialmente adequado para a solução desse problema Rearranjando a Equação 45 para resolvêla para uij 1Δx²ui1j ui1j 1Δy²uij1 uij1 2 1Δx² 1Δy² 413 a qual pode ser utilizada no método iterativo de GaussSiedel Uma estimativa inicial de todos uij é necessária mas isso pode ser facilmente obtido calculando a média das condições de contorno de Dirichlet A convergência do método de GaussSiedel é garantida para um sistema de equações predominantemente diagonal No entanto sua convergência pode ser lenta para a solução de equações de equações diferenciais elípticas O método de sobrerelaxação pode ser utilizado para acelerar a taxa de convergência Essa técnica aplica o seguinte algoritmo ponderador na avaliação dos novos valores de uij em cada iteração do algoritmos de GaussSiedel uijnovo wuijEq 413 1 wuijantigo 414 Devese tomar um cuidado especial no processamento de nós nas fronteiras se uij nesses nós é calculado de uma maneira diferente de diferenças finitas Nessas casos quando é calculado o novo valor de uij a equação apropriada deve ser aplicada em vez da Equação 413 O parâmetro de relaxação pode assumir os valores naos seguinte intervalos 0 w 1 subrelaxação 1 w 2 sobrerelaxação Quando w 1 esse método é exatamente o mesmo que o método GaussSiedel nãomodificado 55 Aula 4 Solução Numérica de Equações Diferenciais Parciais 3 Conteúdo 41 Equações diferenciais parciais elípticas 53 42 Equações diferenciais parciais parabólicas 57 421 Métodos Explícitos 58 422 Métodos Implícitos 61 423 Método das Linhas 62 43 Equações diferenciais Hiperbólicas 64 431 Métodos Explícitos 64 432 Métodos Implícitos 65 44 Fronteiras irregulares e Sistemas com coordenadas polares 66 45 Estabilidade 70 46 Elementos Finitos 72 41 Equações diferenciais parciais elípticas Equações diferenciais parciais elípticas são frequentemente encontradas em operações que envolvem condução de calor e difusão no estado estacionário Por exemplo na condução de calor tridimensional no estado estacionário em sólidos temse a seguinte representação da equação da condução de calor ²Tx² ²Ty² ²Tz² 0 41 De maneira similar a segunda lei de Fick para difusão é simplificada para ²cAx² ²cAy² ²cAz² 0 42 quando é assumido o estado estacionário A discussão das soluções numéricas das equações diferenciais elípticas começa examinandose o problema na forma bidimensional geral Equação de Laplace ²ux² ²uy² 0 43 0 0 0 No caso de uma malha com intervalos equidistantes ou seja quando Δx Δy a Equação 413 é simplificada para uij ui1j ui1j uij1 uij1 4 415 o que demonstra de maneira simples que o valor da variável dependente no ponto pivô ij na equação de Laplace é a média aritmética dos valores nos pontos da malha à esquerda à direita acima e abaixo do ponto pivô Isso é demonstrado pela molécula computacional da figura abaixo a qual é referida por Estrela de 5 pontas A equação diferencial parcial elíptica tridimensional ²ux² ²uy² ²uz² 0 416 pode da mesma maneira ser convertida em equações lineares algébricas utilizando aproximações por diferenças finitas no espaço tridimensional Aplicandose as diferenças finitas na Equação 416 obtémse 1Δx²ui1jk 2uijk ui1jk 1Δy²uij1k 2uijk uij1k 1Δz²uijk1 2uijk uijk1 0 417 Para uma malha equidistante Δx Δy Δz a equação acima se reduz a ui1jk ui1jk uij1k uij1k uijk1 uijk1 6 418 Em paralelo com o caso bidimensional o valor da variável dependente no ponto pivô ijk é a média aritmética dos valores adjacentes do pivô na malha A molécula computacional para a equação diferencial elíptica tridimensional é mostrada na figura abaixo A forma nãohomogênea da Equação de Laplace é a equação de Poisson ²ux² ²uy² fxy 419 e também pertence à classe de equações diferenciais elípticas A forma da equação de Poisson ²Tx² ²Ty² Qxyk 420 é utilizada para descrever a condução de calor em uma placa sólida bidimensional como uma fonte interna de calor Qxy é o calor gerado por unidade de volume por tempo e k é a condutividade térmica do material A formulação da equação de Poisson em diferenças finitas é 21Δx² 1Δy²uij 1Δx²ui1j 1Δx²ui1j 1Δy²uij1 1Δy²uij1 fij 421 ou uij 1Δx²ui1j ui1j 1Δy²uij1 uij1 21Δx² 1Δy² fij 21Δx² 1Δy² 422 42 Equações diferenciais parciais parabólicas Exemplos clássicos de equações diferenciais parabólicas são a equação de condução de calor dinâmica α²Tx² ²Ty² ²Tz² Tt 423 E a segunda lei de Fick para a difusão DAB²cAx² ²cAy² ²cAz² cAt 424 Com condições de contorno de Dririchlet Neumann ou Cauchy Considerase primeiramente essa classe de equação na forma geral unidimensional ut α ²ux² 425 Aqui serão desenvolvidos métodos de solução da equação acima utilizando diferenças finitas 421 Métodos Explícitos As derivadas são expressas em termos das diferenças centrais em torno do ponto in utilizando o contador i para a direção x e o contador n para a direção t ²ux²in 1Δx²ui1n 2uin ui1n OΔx² 426 utin 12Δtuin1 uin1 OΔt² 427 Combinandose as Equações 425 426 e 427 e rearranjando uin1 uin1 2αΔtΔx²ui1n 2uin ui1 OΔx² Δy² 428 Essa é uma fórmula algébrica explícita que calcula o valor da variável dependente no próximo passo de tempo ujn1 a partir dos valores do tempo atual e tempos passados Uma vez que o valor inicial e as condições de contorno do problema estão especificadas a solução da fórmula explícita é direta Contudo essa fórmula explícita em particular é instável pois contém termos negativos no lado direito¹ Como uma regra geral quando todos os valores são arranjados do lado direito na formulação da equação das diferenças finitas se qualquer um dos coeficientes for negativo a solução é instável Essa afirmação é apresentada mais precisamente pela regra da positividade Para a equação uin1 Aui1n Buin Cui1n 429 se A B e C são positivos e A B C 1 então o esquema numérico é estável Com a finalidade de se eliminar o problema de instabilidade a derivada de primeira ordem é substituída na Equação 425 pela diferença progressiva utin 1Δtuin1 uin OΔt 430 Combinando as Equações 425 426 e 430 obtémse a fórmula explícita uin1 αΔtΔx²ui1n 1 2αΔtΔx²uin αΔtΔx²ui1n OΔx² Δt 431 ¹o que será tratado mais adiante no curso Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda Essa é chamada aproximação implícita regressiva que também pode ser obtida pela aproximação de primeira ordem da derivada utilizando diferenças regressivas em in1 e a aproximação de segunda ordem para a derivada pela diferença central em in1 Finalmente quando θ 12 a Equação 448 resulta na fórmula implícita de CrankNicolson αΔtΔx² ui1n1 2 1 αΔtΔx² uin1 αΔtΔx² ui1n1 αΔtΔx² ui1n 2 1 αΔtΔx² uin αΔtΔx² ui1n 450 Para uma solução implícita para a equação parabólica nãohomogênea ut α ²ux² fxt 451 Pelo método acima devese calcular o valor de f no ponto médio in12 a qual é aproximada pela média dos valores de f nos pontos da malha in1 e in fin12 12 fin1 fin 452 Colocandose as Equações 446 447 considerando θ 12 e 452 na Equação 436 resulta em αΔtΔx² ui1n1 2 1 αΔtΔx² uin1 αΔtΔx² ui1n1 Δtfin1 αΔtΔx² ui1n 2 1 αΔtΔx² uin αΔtΔx² ui1n Δtfin 453 A Equação acima é conhecida como a fórmula de CrankNicolson implícita para a solução de equações diferenciais parciais parabólicas nãohomogêneas Quando escrita para toda a malha das diferenças as fórmulas implícitas geram conjuntos de equações algébricas lineares cuja matriz dos coeficientes é geralmente uma matriz tridiagonal Esse tipo de problema pode ser resolvido usando uma eliminação de Gauss ou mais eficientemente utilizando o algoritmo de Thomas variação da eliminação de Gauss Fórmulas implícitas do tipo descrito acima são incondicionalmente estáveis Isso pode ser generalizado para a maioria das aproximações explícitas por diferenças finitas são condicionalmente estáveis enquanto que a maioria das aproximações implícitas são incondicionalmente estáveis Os métodos explícitos no entanto são computacionalmente mais fáceis de resolver do que as técnicas implícitas 423 Método das Linhas Outra técnica Para a solução de equações diferenciais parabólicas é o método das linhas Ele é baseado no conceito do converter a equação diferencial parcial em uma série de equações diferenciais ordinárias discretizando somente as derivadas espaciais usando diferenças finitas e deixando as derivadas em relação ao tempo inalteradas Esse conceito aplicado à Equação 425 resulta em 62 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda duidt αΔx² ui1 2ui ui1 454 Haverá tantas dessas equações diferenciais ordinárias quanto pontos na direção x O conjunto completo dessas equações para 0 i N deve ser du0dt αΔx² u1 2u0 u1 duidt αΔx² ui1 2ui ui1 duNdt αΔx² uN1 2uN uN1 455 As duas equações nas fronteiras tem de ser modificadas de acordo com as condições de contorno especificadas para o problema em particular Por exemplo se uma condição de contorno de Dirichlet é dada em x 0 e t 0 isto é u0 β constante para t 0 456 A primeira equação do conjunto de equações diferenciais ordinárias é modificada para du0dt 0 u00 β 457 Por outro lado se é dada uma condição de Neumann nessa fronteira ou seja ux0t 0 em x 0 e t 0 458 a derivada parcial é substituída por uma aproximação pela diferençacentral ux0t u1 u1 2Δx 0 459 e a primeira equação do conjunto de equações diferenciais ordinárias 455 se torna du0dt αΔx² 2u1 2u0 460 63 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda O conjunto completo de equações diferenciais simultâneas deve ser integrado na direção do avanço do tempo começando com as condições iniciais do problema Esse método resulta em soluções estáveis para equações diferenciais parciais 43 Equações diferenciais Hiperbólicas Equações diferenciais parciais de segunda ordem do tipo hiperbólico ocorrem principalmente em problemas físicos ligados a processos vibracionais Por exemplo a equação unidimensional da onda ρ ²ut² T0 ²ux² fxt 461 descreve o movimento transversal de uma corda vibrando sujeita a uma tensão T0 e uma força externa fxt No caso de uma densidade constante ρ a equação é escrita na forma ²ut² a² ²ux² Fxt 462 onde a² T0ρ Fxt 1ρ fxt Se nenhuma força externa atua sobre a corda a Equação 462 se torna uma equação homogênea ²ut² a² ²ux² 463 A extensão bidimensional da Equação 462 é ²ut² a² ²ux² ²uy² Fxyt 464 a qual descreve a vibração de uma membrana sujeita a um tensão T0 e uma força externa fxyt 431 Métodos Explícitos Para encontrar a solução numérica da equação 463 expandese cada um dos termos das derivadas de segunda ordem em termos das diferenças finitas centrais para obter uin1 uin uin1Δt² a² ui1n 2uin ui1nΔx² OΔx² Δt² 465 Rearranjando para resolver para uin1 uin1 2 1 a²Δt²Δx² uin a²Δt²Δx² ui1n ui1n uin1 OΔx² Δt² 466 64 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda Essa é a solução numérica explícita para a equação hiperbólica 463 A regra da positividade aula passada aplicada à Equação 466 mostra que essa solução é estável se a o seguinte limite imposto pela desigualdade é obedecido a²Δt²Δx² 1 467 De maneira similar a forma homogênea da equação hiperbólica bidimensional ²ut² a²²ux² ²uy² 468 é expandida utilizando as diferenças finitas centrais para resultar em uijn1 2uijn uijn1Δt² a²ui1j 2uijn ui1jnΔx² a²uij1n 2uijn uij1nΔy² OΔx² Δy² Δt² 469 Rearranjando essa equação para a forma explícita usando uma malha equidistante nas direções x e y resulta em uijn1 21 2a²Δt²Δx²uijn uijn1 a²Δt²Δx²ui1jn ui1jn uij1n uij1n 470 Essa solução é estável quando a²Δt²Δx² 12 471 432 Métodos Implícitos Métodos implícitos para a solução de equações diferenciais parciais hiperbólicas podem ser desenvolvidas utilizando a abordagem da variávelponderadora onde o espaço entre as derivadas parciais é ponderado em n1 n e n1 A formulação implícita da Equação 463 é uin1 2uin uin1Δt² a²Δx²θui1n1 2uin1 ui1n1 1 2θui1n 2uin ui1n θui1n1 2uin1 ui1n1 472 onde 0 θ 1 Quando θ 0 Equação 472 é convertida no método explícito 466 Quando θ 12 a Equação 472 tornase uma aproximação do tipo CrankNicolson Métodos implícitos resultam em um conjunto de equações lineares algébricas para as quais a solução pode ser obtida utilizando os métodos de eliminação de Gauss Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda 44 Fronteiras irregulares e Sistemas com coordenadas polares As aproximações por diferenças finitas de equações diferenciais parciais desenvolvidas até o momento são baseadas no sistema de coordenadas cartesianas Com certa frequência no entanto os objetos dos quais as propriedades estão sendo modeladas por equações diferenciais parciais podem ter formas circulares cilíndricas ou esféricas ou fronteiras completamente irregulares As aproximações por diferenças finitas podem ser modificadas para lidar com tais casos Primeiramente considerase um objeto o qual é bem descrito por coordenadas cartesianas em toda sua extensão com exceção das proximidades da fronteira a qual possui formato irregular como mostrado na figura abaixo Há dois métodos para se tratar uma fronteira curvada Um método simples é remodelar a fronteira para passar através do ponto da malha mais próximo a ela Por exemplo na figura acima podese assumir o ponto i j na condição de contorno ao invés do ponto i1 j na fronteira original Apesar de um método simples a aproximação da condição de contorno insere um erro nos cálculos especialmente na fronteira Um método mais preciso de expressar a equação das diferenças finitas em uma condição de contorno irregular é modificála adequadamente Podese usar uma expansão em série de Taylor do valor dependente no ponto i j na direção x para chegar a ui1j uij αΔx uxij α²Δx2 ²ux²ij OΔx³ 473 ui1j uij Δx uxij Δx2 ²ux²ij OΔx³ 474 Eliminando ²ux² das Equações 473 e 474 temse como resultado uxij 1Δx 1α1 αui1j 1 α²uij α²ui1j 475 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda E eliminando ux das Equações 473 e 474 resulta em ²ux²ij 1Δx²2α1 αui1j 1 αuij αui1j 476 De maneira similar na direção y uyij 1Δy 1β1 βuij1 1 β²uij β²uij1 477 ²uy²ij 1Δy²2β1 βuij1 1 βuij βuij1 478 Quando α β 1 as Equações 475478 se tornam idênticas àquelas desenvolvidas anteriormente para o sistema regular de coordenadas cartesianas Portanto para objetos com fronteiras irregulares as equações diferenciais parciais devem ser convertidas a equações algébricas usando Eqs 475478 Para pontos adjacentes à fronteira os parâmetros α e β devem assumir valores que reflitam a forma irregular da fronteira e para os pontos internos distantes da fronteira os valores de α e β devem ser iguais à unidade As Equações 475478 podem ser utilizadas nas fronteiras com condições de contorno de Dirichlet onde a variável dependente é conhecida O tratamento das condições de contorno de Neumann e Robbins onde a derivada normal na fronteira curvada ou irregular é especificada é mais complicado Considerando novamente a figura anterior a derivada normal da variável dependente na fronteira pode ser expressa como uni1j uxi1j cos γ uyi1j sin γ 479 Onde n é o vetor unitário normal à fronteira e γ é o ângulo entre o vetor n e o eixo x As derivadas com respeito a x e y na Equação 479 podem ser aproximadas por expansões em série de Taylor uxi1j uxij αΔx ²ux²ij 480 uyi1j uyij αΔx ²uxyij 481 As derivadas no ponto da malha i j deve ser conhecido com a finalidade de se calcular a derivada normal na fronteira Para a configuração em particular na figura utilizase as diferenças finitas regressivas para avaliar essas derivadas nas Equações 480 e 481 uxij 1Δx uij ui1j 482 ²ux²ij 1Δx²uij 2ui1j ui2j 483 uyij 1Δy uij uij1 484 ²uxyij 1ΔxΔy uij ui1j uij1 ui1j1 485 Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda Combinando as Equações 480 482 e 483 temse que ux i1j 1Δx 1αuij 12αui1j αui2j 486 e a combinação das equações 481 484 e 485 resulta em uy i1j 1Δy 1αuij αui1j 1αuij1 αui1j1 487 Substituindo as Equações 486 e 487 na Equação 479 fornece ao problema a derivada normal a qual pode ser utilizada quando trabalhando com condições de contorno de Neumann ou Robbins em fronteiras irregulares De maneira similar na direção y un ij1 ux ij1 cos γ uy ij1 sin γ 488 onde ux ij1 1Δx 1βuij βuij1 1βui1j βui1j1 489 uy ij1 1Δy 1βuij 12βuij1 αuij2 490 491 É importante lembrar que as Equações 486 487 489 e 490 são específicas para a configuração mostrada na figura para outras possíveis configurações diferenças progressivas ou regressivas ou uma combinação de ambas em diferentes direções podem ser utilizadas para tratar da derivada na condição de contorno Objetos de formato cilíndrico são mais convenientemente representados em forma de coordenadas polares A transformação de uma coordenada cartesiana para o sistema de coordenadas polares é feito utilizando a seguinte representação que é baseada na figura abaixo Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda x r cos θ y r sin θ r x²y² θ tan¹yx 492 O operador laplaciano em coordenadas polares se torna ²ux² ²uy² ²ur² 1r ur 1r² ²uθ² 493 A segunda lei da difusão de Fick em coordenadas polares fica cAt DAB ²cAr² 1r cAr 1r² ²cAθ² ²cAz² 494 Usando a malha das diferenças finitas para as coordenadas polares mostradas na Figura figura as derivadas parciais são aproximadas por ²ur² ij 1Δr² uij1 2uij uij1 495 ²uθ² ij 1Δθ² ui1j 2uij ui1j 496 ur ij 12Δr uij1 uij1 497 onde j e i são contadores nas direções r e θ respectivamente Derivadas parciais nas dimensões z e t não mostradas na figura podem ser expressas de maneira similar através do uso de subscritos adicionais Equações diferenciais parciais nãolineares A discussão desenvolvida nesse curso foi focada em equações diferenciais parciais que resultam em conjuntos de equações algébricas lineares quando expressas em aproximações por diferenças finitas Por outro lado se a equação diferencial parcial é nãolinear por exemplo Modelagem e Simulação de Processos II Prof Dr Júlio Cesar de Carvalho Miranda u ²ux² u ²uy² fu 498 A discretização resultante da aproximação por diferenças finitas pode gerar conjuntos de equações algébricas nãolineares A solução desse tipo de problema pode demandar a aplicação do método de Newton curso de Modelagem I para equações nãolineares simultâneas 45 Estabilidade Nessa aula serão discutidos alguns aspectos da estabilidade das aproximações por diferenças finitas utilizando o conhecido procedimento de von Neumann Esse método introduz um erro inicial representado por uma série de Fourier finita e verifica como esse erro se propaga durante a solução O método de von Neumann se aplica a problemas de valor inicial por essa razão é utilizado para analisar a estabilidade do método explícito para equações parabólicas e também o método explícito para equações hiperbólicas desenvolvidos anteriormente Definese o erro ϵmn como a diferença entre a solução umn da aproximação por diferença finita e a solução exata ūmn da equação diferencial no passo mn ϵmn umn ūmn 499 A solução por diferenças finitas explícitas da equação diferencial parabólica pode se escrita para umn1 e ūmn1 como a seguir umn1 α Δt Δx² um1n 1 2 α Δt Δx² umn α Δt Δx² um1n REmn1 4100 ūmn1 α Δt Δx² ūm1n 1 2 α Δt Δx² ūmn α Δt Δx² ūm1n TEmn1 4101 4102 onde REmn1 e TEmn1 são respectivamente os erros de arredondamento no passo mn1 Combinando as Equações 4994101 obtémse ϵmn1 α ΔtΔx² ϵm1n 1 2 α ΔtΔx² ϵmn α ΔtΔx² ϵm1n REmn1 TEmn1 4103 Essa é uma equação de diferenças finitas nãohomogênea em duas dimensões representando a propagação do erro durante a solução numérica da equação parcial diferencial parabólica A Solução dessa equação de diferenças finitas é bastante difícil de se obter Por essa razão a análise de von Neumann considera somente a parte homogênea da Equação 4103 ϵmn1 α ΔtΔx² ϵm1n 1 2 α ΔtΔx² ϵmn α ΔtΔx² ϵm1n 0 4104 ut cu au f 4125 e ²ut² cu au f 4126 onde d c a e f são funções complexas no domínio da solução e também podem ser funções do tempo O símbolo é o operador diferencial vetorial não confundir com o vetor da diferença regressiva Nas toolboxes as equações 41244126 são descritas como elíptica parabólica e hiperbólica respectivamente independentemente dos valores dos coeficientes e das condições de contorno Bibliografia 1 A Constantinides and N Mostoufi Numerical Methods for Chemical Engineers with MATLAB Applications Prentice Hall PTR 1999 2 S C Chapra and R P Canale Numerical Methods for Engineers McGrawHill 2015