Esquemas Essencialmente Não-Oscilatórios de Alta Precisão ... · Matemática Aplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro (UFRJ), como parte dos

Universidade Federal do Rio de Janeiro

Esquemas Essencialmente Não-Oscilatórios deAlta Precisão para Leis de Conservação

Hiperbólicas

Marcos Castro

Rio de Janeiro

Novembro de 2009


Esquemas Essencialmente Não-Oscilatórios de Alta Precisãopara Leis de Conservação Hiperbólicas

Marcos Castro C. T. de Azevedo

Dissertação de Mestrado apresentada ao Programa de Pós-graduação em MatemáticaAplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro (UFRJ),como parte dos requisitos necessários obtenção do título de Mestre em Matemática Apli-cada.

Orientador: Prof. Bruno Alexandre Soares da Costa.

Rio de JaneiroNovembro de 2009

ii


Esquemas Essencialmente Não-Oscilatórios de Alta Precisãopara Leis de Conservação Hiperbólicas

Marcos Castro C. T. de Azevedo

Dissertação de Mestrado apresentada ao Programa de Pós-graduação em MatemáticaAplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro (UFRJ),como parte dos requisitos necessários obtenção do título de Mestre em Matemática Apli-cada.

Aprovada por:

———————————————————Presidente, Prof. Bruno Alexandre Soares da Costa

———————————————————Prof. Wai Sun Don

———————————————————Prof. Marco Aurélio Palumbo Cabral


iii

Castro, Marcos.Esquemas Essencialmente Não-Oscilatórios de

Alta Precisão para Leis de Conservação Hiperbólicas / MarcosCastro. – Rio de Janeiro: UFRJ/ IM, 2009.

viii, 132f.:il.; 29,7 cm.Orientador: Bruno Alexandre Soares da Costa.Dissertação (Mestrado) — UFRJ/ IM/ Programa de

Pós-graduação em Matemática Aplicada, 2009.Referências Bibliográficas: f. 131-132 .1. Método WENO. 2. Leis deConservação. 3. Métodos Numéricos

I. Costa, Bruno Alexandre Soares.II. Universidade Federal do Rio de Janeiro,Instituto de Matemática. III. Esquemas EssencialmenteNão-Oscilatórios de Alta Precisão para Leisde Conservação Hiperbólicas.

iv

Dedicado à minha famlia.

v

Agradecimentos

Aos meus pais e família, aqueles que tornaram tudo isto possível.

A Luciana, por todo o carinho e apoio, pelos quais não tenho como agradecer atravésde palavras.

Ao meu orientador, Bruno Costa, pelos ensinamentos, paciência e compreensão duranteeste trabalho e os anos que o precederam.

Ao professor Marco Cabral, um dos meus primeiros professores e um dos grandes res-ponsáveis para que eu seguisse este caminho.

Ao professor Wai Sun Don, que, assim como o professor Marco Cabral, prontamenteaceitou o convite para compor a banca avaliadora.

Aos meus amigos e colegas, que contribuíram de diversas maneiras para que eu chegasseaté aqui. Não os citarei nominalmente para não cometer injustiças e não me alongar.

Marcos CastroNovembro de 2009

vi

“Não existe um caminho para a

paz. A paz é o caminho."

Mahatma GandhiLíder espiritual e pacifista indiano

vii

Resumo

Esquemas Essencialmente Não-Oscilatórios de Alta Precisão para Leis de

Conservação Hiperbólicas

Marcos Castro

Resumo da dissertação de Mestrado submetida ao Programa de Pós-graduação emMatemática Aplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro(UFRJ), como parte dos requisitos necessrios obtenção do título de Mestre em MatemáticaAplicada.

Resumo: Nesta dissertação desenvolvemos um novo esquema essencialmentenão-oscilatório ponderado (WENO) de quinta ordem adicionando um indica-dor de suavidade de alta ordem, obtido por uma simples combinação linear(e de baixo custo computacional) dos indicadores de suavidade já existentes.Este novo esquema, denominado WENO-Z, tem um custo computacional equi-valente ao clássico WENO-JS e consideravelmente menor do que o mapeado(WENO-M), uma vez que não há nenhum mapeamento envolvido. Tambémobservamos mais de perto as expansões de Taylor dos polinômios de Lagrangedos subestênceis WENO e as simetrias herdadas dos indicadores de suavidadede ordem menor para obter uma fórmula geral para os indicadores de suavi-dade de alta ordem, permitindo a extensão do método WENO-Z para todas asordens (ímpares) de aproximação. Observamos também a boa aproximação doWENO-Z em pontos críticos de soluções suaves tal como suas característicasnuméricas diferenciadas, derivadas do novo conjunto de pesos não-linares etambém mostramos que em relação à dissipação numérica weno-Z ocupa umaposição intermediária entre WENO-JS e WENO-M.

Palavras–chave., WENO, Leis de Conservação, Indicadores de Suavidade, Pesos Não-Lineares.


viii

Abstract

High Order Weighted Essentially Non-Oscillatory Schemes for Hyperbolic

Conservation Laws

Marcos Castro

Abstract da dissertação de Mestrado submetida ao Programa de Pós-graduao em Ma-temática Aplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro(UFRJ), como parte dos requisitos necessários obtenção do título de Mestre em Matemá-tica Aplicada.

Abstract: In this work we design a new fifth order WENO finite-differencescheme by adding a higher order smoothness indicator which is obtained asa simple and inexpensive linear combination of the already existing low ordersmoothness indicators. Moreover, this new scheme, dubbed as WENO-Z, has aCPU cost which is equivalent to the one of the classical WENO-JS and signifi-cantly lower than that of the mapped WENO-M„ since it involves no mappingof the nonlinear weights. We also take a closer look at Taylor expansions ofthe Lagrangian polynomials of the WENO substencils and the related inhe-rited symmetries of the classical lower order smoothness indicators to obtaina general formula for the higher order smoothness indicators that allows theextension of the WENO-Z scheme to all (odd) orders of accuracy. We furtherinvestigate the improved accuracy of the WENO-Z schemes at critical pointsof smooth solutions as well as their distinct numerical features as a result ofthe new sets of nonlinear weights and we show that regarding the numericaldissipation WENO-Z occupies an intermediary position between WENO-JSand WENO-M.

Keywords.WENO, Conservation Laws, Smoothness Indicators, Non-linear Weights.

Rio de JaneiroNovembro 2009

Sumário

1 Introdução 3

2 Tópicos introdutórios 11

2.1 Aproximação por polinômios . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.1.1 Interpolação via polinômios de Lagrange . . . . . . . . . . . . . . . 14

2.2 Indicadores de suavidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

2.2.1 Um exemplo: a medida da variação total . . . . . . . . . . . . . . . 19

2.3 Leis de conservação hiperbólicas . . . . . . . . . . . . . . . . . . . . . . . . 21

2.3.1 Equações de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4 Métodos numéricos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.4.1 Resolução numérica de leis de conservação . . . . . . . . . . . . . . 29

2.4.2 Aproximação por estêncil fixo . . . . . . . . . . . . . . . . . . . . . 30

2.4.3 O que torna um método conservativo? . . . . . . . . . . . . . . . . 34

2.4.4 Oscilações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3 Métodos essencialmente não-oscilatórios 39

3.1 O método ENO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.2 Esquemas WENO . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.2.1 Pesos ideais e comparação com upstream central . . . . . . . . . . . 47

1

SUMÁRIO SUMÁRIO

3.2.2 WENO clássico: os pesos de Jiang-Shu . . . . . . . . . . . . . . . . 49

3.2.3 Pontos críticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

3.2.4 WENO mapeado: corrigindo algumas falhas . . . . . . . . . . . . . 60

4 O método WENO-Z 65

4.1 Uma nova fórmula para os pesos do WENO . . . . . . . . . . . . . . . . . 66

4.1.1 Generalização para todo r . . . . . . . . . . . . . . . . . . . . . . . 69

4.2 Resultados auxiliares . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.2.1 Representação dos polinômios fk . . . . . . . . . . . . . . . . . . . 71

4.2.2 Ordem dos indicadores de suavidades βk . . . . . . . . . . . . . . . 82

4.3 Existência de τn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

4.4 Ordem máxima de τn . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

4.5 Uma alternativa para obter ordem n . . . . . . . . . . . . . . . . . . . . . 106

5 Resultados numéricos 109

5.1 Análise dos pesos WENO . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.1.1 Ordens mais altas . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

5.2 Tempo de execução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

5.3 Equações de Euler . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

5.3.1 Comparação entre os esquemas WENO . . . . . . . . . . . . . . . . 122

5.3.2 Conclusão . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

2

Capítulo 1

Introdução

Uma equação de lei de conservação hiperbólica é uma equação da forma

∂u

∂t+∇ · F(u) = 0,

onde F é definido como o fluxo conservado.

Sem perda de generalidade, restringiremos nossa discussão ao caso unidimensional,

∂

∂tu(x, t) +

∂

∂xf(u(x, t)) = 0. (1.1)

O nome lei de conservação se deve ao seguinte fato: considere u(x, t) a densidade de uma

certa grandeza, de modo que∫ x2

x1

u(x, t)dx

represente a quantidade desta grandeza no intervalo [x1, x2] no tempo t. Integrando (1.1)

em [x1, x2]× [t1, t2] chegamos à fórmula conservativa

∫ x2

x1

u(x, t2)dx =

∫ x2

x1

u(x, t1)dx+

∫ t2

t1

f(u(x1, t))dx−

∫ t2

t1

f(u(x2, t))dx

3

CAPÍTULO 1. INTRODUÇÃO

Informalmente, a quantidade da grandeza no instante t2 > t1 é igual à quantidade no

instante t1 mais o que entra pela extremidade x1 e menos o que sai pela extremidade x2.

Por isto, a função f recebe o nome de fluxo.

As equações de leis de conservação frequentemente aparecem na física quando lidamos

com quantidades que se conservam. Algumas das equações mais famosas são a equação

do transporte

ut + αux = 0,

a equação de Burgers

ut +

(

u2

2

)

x

= ut + uux = 0

e o sistema de equações de Euler

ρ

ρv

E

t

+

ρv

ρv2 + p

v(E + p)

x

= 0

As equações de Euler têm aplicação na dinâmica de fluidos, especialmente sobre fluidos

em alta velocidade. A variável ρ representa a densidade de massa, v é a velocidade, E

é a energia total e p é a pressão, geralmente uma função das três primeiras variáveis. A

primeira linha do sistema corresponde à conservação de massa, a segunda à conservação

de momento e a terceira à conservação de energia.

As equações de leis de conservação não-lineares em geral não possuem uma solução

analítica explícita. Por este motivo métodos numéricos são frequentemente utilizados para

obter soluções. Porém, pela natureza das leis de conservação, alguns cuidados têm que ser

tomados. Uma equação de lei de conservação pode ter uma solução u(x, t) descontínua

mesmo quando a condição inicial u0(x) é C∞. Estas soluções descontínuas, mais apropri-

4


adamente chamadas de soluções fracas, em geral não são únicas para um dado problema,

apesar de haver apenas uma solução fisicamente relevante (aquela que satisfaz condições

de entropia, adicionais ao problema). Isto torna mais complicada a tarefa de achar uma

solução numérica.

Primeiramente, para obtermos a solução fraca fisicamente relevante o método numérico

deve ser conservativo e as soluções por ele obtidas devem satisfazer a condições de entropia.

Não entraremos em detalhes sobre este importante tópico. Uma boa referência se encontra

em [6, cap. 12].

Em segundo lugar, o método numérico que vamos usar deve se comportar bem diante

de descontinuidades. Em geral, um método comum baseado em volumes ou diferenças

finitas é dissipativo (suavizando as descontinuidades) se sua ordem é 1 ou dispersivo

(com oscilações parecidas com o fenômeno de Gibbs) se a ordem é maior ou igual a 2.

Estas oscilações tornam estes métodos instáveis e as soluções são qualitativamente ruins.

Entretanto, as dissipações geradas por métodos de ordem 1 são severas, de maneira que

muitos pontos são necessários para se obter uma solução confiável.

Se quisermos soluções de ordem alta (de 2 em diante) então precisamos usar métodos

especiais para não gerarmos oscilações. Uma classe de métodos que não são dispersivos

mas sim dissipativos mesmo em ordens altas é a dos chamados métodos Essencialmente

Não-Oscilatórios. O nome se deve ao método denominado essentially non-oscillatory

schemes (conhecido também por sua abreviatura ENO), proposto no artigo [1] de Harten,

Enquist, Osher e Chakravarthy. Mais tarde, outros métodos baseados no ENO surgiram.

O mais importante deles foi o esquema essencialmente não-oscilatório ponderado (WENO,

do original em inglês weighted essentially non-oscillatory), desenvolvido nos artigos [2] de

Liu, Osher e Chan e [4] de Jiang e Shu. Recentemente, o artigo [3] de Henrick, Aslam e

Powers lançou uma nova luz sobre o WENO, com o desenvolvimento do WENO mapeado.

Há ainda outras modificações, como por exemplo oo ENO geométrico [8] e o power ENO

5


Figura 1.1: Solução numérica da equação do transporte obtida com um método de estêncilfixo de ordem 1. Note como o método é dissipativo

[7]. O estudo dos esquemas essencialmente não-oscilatórios, em especial o ENO e os

esquemas baseados no WENO, será o assunto principal deste texto.

O esquema ENO se baseia numa escolha: para cada ponto, há r conjuntos de r pontos

adjacentes (chamados de estênceis) que são candidatos a serem usados para aproximar

f ′. Usando medidores de suavidade, o algoritmo escolhe qual dentre os estênceis é o mais

suave e o utiliza para aproximar f ′, evitando as oscilações.

Entretanto, em regiões onde a solução é suave tem um inconveniente: desprezamos

informação útil sobre a função quando usamos 1 estêncil na aproximação e descartamos os

outros r− 1. Se usássemos todos os estênceis, teríamos uma quantidade maior de pontos

(2r − 1, para ser mais preciso) disponíveis para a aproximação e, consequentemente,

aumentaríamos a ordem de convergência do método. O esquema WENO faz isto ao

atribuir pesos a cada estêncil. A aproximação final é obtida pela soma ponderada da

aproximação em cada um dos estênceis. Os pesos foram desenvolvidos de tal forma que

um estêncil onde a função é descontínua recebe peso praticamente nulo. Deste modo, o

6


Figura 1.2: Solução numérica da equação do transporte obtida com um método de estêncilfixo de ordem 5. As oscilações vistas são comuns em métodos alta ordem

WENO emula o comportamento do ENO nas vizinhanças de descontinuidades, evitando

assim oscilações na solução. Ao mesmo tempo, em regiões contínuas todos os estênceis

são levados em conta, permitindo um ganho em termos de ordem de convergência. Os

métodos essencialmente não-oscilatórios se encontram em uma página recente da história

da matemática (o artigo pioneiro [1]desenvolvido neste campo data de 1987).

No capítulo 2 apresentaremos tópicos preliminares, necessários para o entendimento

dos métodos ENO. Na primeira seção deste capítulo estudaremos a aproximação polino-

mial, através da interpolação de Lagrange. Os métodos de diferença finita que veremos,

incluindo os essencialmente não-oscilatórios, utilizam este tipo de aproximação polinomial

em suas fórmulas.

Na segunda seção do capítulo 2 entraremos em detalhes sobre os medidores de suavi-

7


dade - outra peça fundamental para os esquemas ENO. Definiremos o que é um medidor

de suavidade; informalmente falando, trata-se de um operador que associa um valor a um

conjunto de valores pontuais de uma função. A magnitude deste valor indica o grau de

suavidade de uma função: quanto menor for o valor, mais suave será a função no domínio

da medição. Existem diversos tipos de medidores de suavidade, mas estudaremos apenas

a medida da variação total, que será o medidor utilizado nos esquemas ENO.

Na terceira seção do capítulo 2 abordaremos uma importante propriedade das leis de

conservação hiperbólicas, que é a possibilidade de diagonalizar o sistema de equações, de

modo a desacoplar as variáveis, isto é, fazer com que cada equação dependa de apenas

uma variável. Veremos que esta propriedade é fundamental para a resolução de certas

equações, como o sistema de equações de Euler, por exemplo.

Na quarta seção do capítulo 2 veremos questões relevantes a esquemas numéricos

de resolução de leis de conservação, como as condições para estabilidade e a propriedade

conservativa de um método. Nesta seção, veremos ainda os chamados esquemas de estêncil

fixo e mostraremos seu comportamento oscilatório. Os esquemas de estêncil fixo recebem

esta terminologia devido ao fato de que o estêncil utilizado não varia em função do ponto

da malha onde se faz a aproximação. Estes esquemas são bem simples de se implementar;

porém, eles geram oscilações significativas nas proximidades de uma descontinuidade.

Apesar deste mal comportamento diante deste tipo de problema, o estudo dos esquemas

de estêncil fixo se faz necessário, já que o ENO pode ser entendido como uma modificação

deste tipo de esquema.

No capítulo 3 estudaremos alguns tipos de métodos ENO: o ENO original [1], o WENO

clássico com os pesos de Jiang-Shu [4] e o WENO mapeado, proposto em [3] para contornar

problemas encontrados em vizinhanças de pontos críticos (isto é, pontos de derivada zero).

O nome WENO mapeado se deve ao fato de que um mapeamento é aplicado aos pesos de

Jiang-Shu de modo a corrigir estes problemas observados. Além disso, este mapeamento

8


faz com que os pesos calculados perto das descontinuidades ganhem mais importância,

tornando o método menos dissipativo. No entanto, vemos que o mapeamento acaba por

deixar o esquema mais caro computacionalmente.

No capítulo 4 introduziremos o método WENO-Z [13], que introduz um novo cálculo

para os pesos de Jiang-Shu: a ideia é inserir na formulação dos pesos um medidor de

suavidade global que contém informações sobre todos os r estênceis. Observamos que o

WENO-Z sempre alcança ordem ótima em regiões contínuas afastadas de pontos críticos

e que, em ordens mais altas, a optimalidade ainda é alcançada mesmo com a presença

destes pontos críticos.

Finalmente, no capítulo 5 mostraremos soluções numéricas de alguns exemplos clássi-

cos, como a equação do transporte e o sistema de equações de Euler, comparando os três

tipos de WENO vistos neste texto. Os resultados encontrados confirmaram que o WENO

mapeado é o menos dissipativo dentre os três esquemas, seguido do WENO-Z e por fim

o WENO clássico. Em termos de desempenho computacional, o WENO-Z, por sua vez,

apresentou vantagens em relação ao WENO mapeado.

9


10

Capítulo 2

Tópicos introdutórios

Antes de iniciarmos o estudo sobre esquemas ENO será preciso examinar dois assuntos

importantes, ainda que relativamente simples. O primeiro deles surge no processo de

aproximação dos esquemas de diferença finita. Neste processo, em um dado momento

estamos interessados em obter coeficientes c−l, c−l+1, . . . , ck que servirão para obtermos

uma aproximação para f ′(xi),, de ordem r, usando os valores de f em um conjunto de r

pontos que inclui xi,

S = xi−l, xi−l+1, . . . , xi, . . . , xi+k , r = k + l + 1,

onde

xi+j = xi + j∆x.

A aproximação assume a forma

∑k

j=−l cjf(xi+j)−∑k

j=−l cjf(xi+j−1)

∆x= f ′(xi) +O(∆xr). (2.1)

11

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS

Há uma formulação equivalente para este problema, na qual é mais simples encontrar

estes coeficientes. Considere o fluxo numérico h definido implicitamente como

f(x) =1

∆x

∫ x+∆x2

x−∆x2

h(ξ)dξ. (2.2)

Esta função satisfaz∂f

∂x

∣

∣

∣

∣

x=xi

=hi+ 1

2− hi− 1

2

∆x. (2.3)

Pode-se mostrar (vide [5]) que os coeficientes que satisfazem

k∑

j=−l

cjf(xi+j) = h(xi+ 12) +O(∆xr) (2.4)

são os mesmos que satisfazem (2.1). Destacamos que para cada j o valor f(xi+j) é

uma média celular de h no intervalo [xi+j− 12, xi+j+ 1

2]. Portanto, o problema que estamos

resolvendo é

Problema 1. Dado os valores f(xi+j) das médias celulares de h, ache coeficientes c−l, c−l+1, . . . , ck

que satisfazem (2.4).

Este problema é facilmente resolvido através da interpolação polinomial. Veremos a

solução na seção seguinte, via polinômios de Lagrange.

O segundo assunto deste capítulo se refere aos chamados medidores de suavidade. Ba-

sicamente, um medidor de suavidade é um operador capaz de determinar a suavidade de

uma função f em um dado intervalo [x0, xr] a partir de um conjunto de valores pontuais

f(x0), f(x1), . . . , f(xr), xi ∈ [x0, xr]. Por causa desta propriedade, os medidores de suavi-

dade são empregados nos esquemas baseados no ENO. Estes esquemas determinam uma

região onde a função é contínua e a utilizam para realizar a aproximação de f ′(xi). Dentre

12

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS2.1. APROXIMAÇÃO POR POLINÔMIOS

diversos exemplos de medidores de suavidade, veremos a medida de variação total, um

medidor introduzido em [4] que é bastante utilizado nos esquemas WENO.

Finalmente, na terceira seção deste capítulo estudaremos os rudimentos da resolução

numérica de equações de leis de conservação. Veremos as condições que o método ne-

cessariamente precisa satisfazer para que tenhamos uma solução correta: a condição de

Courant-Friedrichs-Lax (CFL). Veremos também um método numérico de resolução de

leis de conservação bastante simples. Este esquema, ainda que oscilatório, constitui uma

base da qual os métodos essencialmente não-oscilatórios são uma modificação.

2.1 Aproximação por polinômios

Considere uma função f e uma malha de pontos

S = x0, . . . , xk , x0 < x1 < . . . < xk. (2.5)

Vamos supor que saibamos o valor de f nos pontos da malha

f(xi) = fi, i = 0, . . . , k

mas que não saibamos, ou não queiramos calcular, o valor de f em nenhum outro ponto.

Porém, estamos interessados em obter um valor aproximado de f(x) para x ∈ (x0, xk).

Uma maneira muito natural de se fazer isto é achar o polinômio p(x) que interpola f na

malha S, isto é, o polinômio que satisfaz

p(xi) = fi, i = 0, . . . , k (2.6)

A partir daí, usamos p(x) como uma aproximação para f(x) no intervalor (x0, xk).

13

2.1. APROXIMAÇÃO POR POLINÔMIOSCAPÍTULO 2. TÓPICOS INTRODUTÓRIOS

A escolha do grau de p é muito importante. Se o grau for pequeno pode ser impossível

achar p que satisfaça (2.6). Por outro lado, se o grau for grande demais então pode existir

mais de um polinômio que satisfaça (2.6), e neste caso a aproximação fica igualmente

indeterminada. Felizmente há um fato que nos ajudará a resolver esta questão (omitiremos

a demonstração, apesar de ser relativamente simples).

Teorema 1. Considere uma malha de pontos e uma função f definidas como acima.

Existe um único polinômio p(x) de grau menor ou igual a k tal que p(xi) = f(xi) para

i = 0, ..., k.

Algo também importante é a ordem do erro que obtemos com a aproximação. O

teorema a seguir nos dá esta informação:

Teorema 2. Seja p o único polinômio de grau menor ou igual a k que interpola f na

malha S. Além disto, vamos definir

∆x ≡ max1≤i≤k

xi − xi−1. (2.7)

Se f ∈ Ck[x0, xk] então

f(x)− p(x) = O(∆xk+1)

para todo x ∈ [x0, xk].

2.1.1 Interpolação via polinômios de Lagrange

Há uma maneira bastante simples de calcular o polinômio interpolante p(x). Considere

a família de polinômios li(x) de grau k que é definida como

li(x) =k∏

j=0j 6=i

x− xj

xi − xj

, i = 0, . . . , k

14

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS2.1. APROXIMAÇÃO POR POLINÔMIOS

Os polinômios li(x) são bem definidos desde que a malha só contenha pontos distintos,

o que é o caso de (2.5).

Observe que, por definição,

li(xj) =

0, se i 6= j

1, se i = j(2.8)

Daí segue um fato interessante a respeito destes polinómios:

Teorema 3.∑k

i=0 li(x) = 1 para todo x ∈ R.

Demonstração. Vamos chamar q(x) =∑k

i=0 li(x). Usando (2.8) podemos ver que

q(xj) =k∑

i=0

li(xj) = 1, j = 0, . . . , k.

Portanto, o polinomio r(x) = 1 interpola q(x) nos pontos x0, . . . , xk. Por outro lado, q(x) é

um polinômio de grau no máximo k e, obviamente, ele também se interpola em x0, . . . , xk.

Pelo Teorema 1, existe um único polinômio de grau no máximo k que interpola q(x) em

x0, . . . , xk. Logo∑k

i=0 li(x) ≡ r(x) ≡ 1.

Vejamos como podemos usar os polinomios li(x) para acharmos uma fórmula para o

polinomio interpolante p(x). Usando (2.8), verifica-se facilmente que∑k

i=0 fili(xj) = fj.

Daí segue imediatamente que

p(x) =k∑

i=0

fili(x) =k∑

i=0

fi

k∏

j=0j 6=i

x− xj

xi − xj

(2.9)

é o polinomio que procurávamos. Note que p(x) é uma soma de polinomios de ordem k

e que, por isto, o grau de p(x) é menor ou igual a k. Assim, acabamos de demonstrar a

parte da existência do polinomio interpolante do teorema 1.

15

2.2. INDICADORES DE SUAVIDADECAPÍTULO 2. TÓPICOS INTRODUTÓRIOS

Esta maneira de escrever o polinomio interpolante é conhecida como fórmula interpola-

dora de Lagrange. Há outras maneiras de se calcular o polinômio interpolante, como pela

fórmula de Newton (não abordaremos este assunto aqui). Porém a fórmula de Lagrange

será mais conveniente para nós neste momento. Vejamos o porquê.

Neste texto, em vários momentos, estaremos interessados em obter uma expressão

para a aproximação polinomial de f(x) em termos dos valores conhecidos f0, . . . , fk. Em

outras palavras, vamos querer achar os coeficientes c0, . . . , ck que satisfaçam

f(x) ≈ p(x) =k∑

i=0

cifi (2.10)

Podemos ver, por (2.9), que no contexto atual

ci = li(x) =k∏

j=0j 6=i

x− xj

xi − xj

.

2.2 Indicadores de suavidade

Como saber se uma função f é contínua? Esta pergunta, aparentemente inocente, es-

conde um problema dos mais relevantes. Claro, quando sabemos a definição da função(f(x) =

sin(x) ou f(x) = 1/ (x+ 3), por exemplo), a resposta é mais fácil. Mas e se não tivermos

tanta informação sobre a função? Na prática, quando lidamos com funções num compu-

tador raramente sabemos a definição delas, a não ser que estejamos usando um programa

de matemática simbólica, como o Maple. Por exemplo, em geral não sabemos nada sobre

a solução de uma EDP. Ao invés disto, quando a resolvemos numericamente o que temos

são apenas avaliações pontuais da solução u (uma aproximação dela, na verdade) em um

conjunto finito e discreto de pontos x0, . . . , xk. Será que conhecendo apenas estes va-

lores pontuais é possível afirmar algo a respeito da continuidade de u? Ou, formulando

16

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS2.2. INDICADORES DE SUAVIDADE

melhor,é possível dizer se uma função f(x) é contínua conhecendo-se apenas um conjunto

discreto de valores f(x0), f(x1), . . . , f(xk))? Ou simplesmente, é possível dizer se uma

função f(x) é contínua apenas olhando para o seu gráfico? Qualquer um responderia

ao problema acima dizendo que sim, é possível dizer se uma função é contínua ou não

só olhando para seu gráfico. Isto porque o cérebro humano é equipado com um “sensor

de continuidade” que faz toda a tarefa de detectar descontinuidades automaticamente.

Sabemos quando uma função é suave e quando ela dá saltos. Um olho mais treinado é

capaz até de detectar descontinuidades na segunda derivada de uma função só de olhar

para o seu gráfico.

Mas acontece que no mundo real em que vivemos é impossível representar exatamente

o gráfico de uma função cujo domínio é um conjunto denso. Por exemplo, nunca consegui-

remos desenhar perfeitamente uma reta densa, aquela em que entre dois pontos quaisquer

sempre existe mais um ponto. Isto porque quando desenhamos um gráfico com um lápis,

a linha que parece cheia é na verdade uma cadeia discreta de moléculas de grafite. Ou

então quando desenhamos o gráfico no computador, usamos um conjunto finito de pixels,

e assim por diante. Porém, quando olhamos para um gráfico, nosso cérebro não é capaz

de detectar os buracos entre uma molécula e outra (ou entre um pixel e outro), criando

a ilusão de que vemos de fato uma reta densa. E por causa disto que o nosso “sensor de

continuidade” ainda funciona mesmo nestes gráficos que não representam tão exatamente

a função.

Observe que o gráfico de f(x) na verdade ilustra a avaliação de f em um conjunto finito

e discreto de pontos. Por este motivo, os dois problemas acima são semelhantes. Logo,

já que a nossa mente consegue olhar para um gráfico esburacado e dizer se uma função é

contínua, então pode haver alguma maneira de se fazer o mesmo matematicamente. E há

uma maneira, de fato: usando-se os medidores de suavidade, que definimos a seguir:

17


Definição 1. Considere uma malha de pontos x0 < x1 < . . . < xk tal que ∆x =

max1≤i≤k xi − xi−1 < 1. Além disto, considere uma função f : D ⊂ R −→ R con-

tínua por partes em (x0, xk) tal que, para todo x ∈ D, temos f(x) = limt→x+ f(t) ou

f(x) = limt→x− f(t). Seja f(xi) = fi, i = 0, . . . , k. Um medidor de suavidade β é um

operador que, a partir dos valores f0, f1, . . . , fk, é capaz de avaliar a suavidade de f em

(x0, xk) no seguinte sentido:

• Se f é contínua em (x0, xk) então β(f0, f1, . . . , fk) = O(∆xp) para algum p ≥ 1;

• Se f é descontínua em (x0, xk) então β(f0, f1, . . . , fk) = O( 1∆xq ) para algum q ≥ 0;

O nome "medidor de suavidade"vem do fato de que quanto mais suave for f em (x0, xk)

em geral menores serão os valores de β(f0, f1, . . . , fk). Deste modo é possível detectar a

presença de descontinuidades em f num domínio apenas olhando para a magnitude de β.

Podemos ver que nesta definição é exigido que f seja contínua por partes. Vejamos o

porquê. Para uma função f descontínua em um conjunto de medida não-nula pode ser

impossível desenharmos seu gráfico. Seria preciso que utilizássemos um instrumento com

resolução infinita, capaz de realizar a impossível tarefa de mostrar o valor da função em

todos os pontos onde há descontinuidade. Do mesmo modo, seria preciso um ∆x infini-

tesimalmente pequeno para que pudéssemos captar todas as nuances de f neste conjunto

de medida não-nula, o que é impossível neste universo em que vivemos. Além do mais,

nas aplicações com maior interesse prático é muito incomum lidarmos com funções que

não sejam contínuas por partes, fato que torna esta imposição algo nem tanto restritivo.

Mas ainda há outro problema. Considere a função

g(x) =

1, se x 6= 0

0, se x = 0(2.11)

18

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS2.2. INDICADORES DE SUAVIDADE

Neste caso, qualquer medidor de suavidade definido numa malha que não contenha o

ponto x = 0 não será capaz de detectar a descontinuidade de g(x), apesar de g(x) ser

contínua por partes. Por causa de casos assim é que exigimos que f(x) = limt→x+ f(t) ou

f(x) = limt→x− f(t) para todo x ∈ D.

A definição (1) é genérica o suficiente para englobar várias classes de medidores de

suavidade. De fato, há medidores de suavidade baseados na interpolação polinomial,

como a comparação com a aproximação polinomial e as diferenças divididas de Newton;

há a medida da variação total, que, como o nome diz, usa a variação total da função; há

também os medidores que se utilizam das séries de Fourier da função [10], [11] e medidores

de suavidade que se baseiam em outros artifícios matemáticos (por exemplo, [12]). Por

brevidade, apenas abordaremos a medida de variação total.

2.2.1 Um exemplo: a medida da variação total

Introduzida por Jiang e Shu [4], a medida de variação total é o medidor de suavidade

padrão para o WENO. Ela é baseada no conceito de variação total, ou TV. A variação

total de uma função f no intervalo [a, b], é definida como

TV [f ] = supP

∑

i|xi,xi+1∈P

|f(xi+1)− f(xi)| ,

onde P = a = x0 < x1 < . . . < xN = b é uma partição de [a, b]. Pode-se mostrar que,

se f é contínua por partes em [a, b], então

TV [f ] =

∫ b

a

|f ′(x)| dx = ‖f ′(x)‖L1 .

19


Aqui vamos usar a versão L2 da variação total,

TVL2 [f ] =

∫ b

a

(f ′(x))2dx = ‖f ′(x)‖

2L2 .

Considere o estêncil

Sk = xi+k−r+1, xi+k−r+2, . . . , xi+k , k = 0, 1, . . . , r − 1

e seja pk(x) o polinômio (de ordem até r−1) que interpola f em Sk. A medida de variação

total de ordem r de f em Sk é definida como

βk =r−1∑

n=1

∆x2n−1TVL2 [p(n)k ]

=r−1∑

n=1

∫ xi+1

2

xi− 1

2

∆x2n−1(

p(n)k (x)

)2

dx,

onde p(n)k denota a n-ésima derivada de pk(x). Ou seja, βk é em essência uma aproximação

da soma das variações totais das derivadas de f . O termo ∆x2n−1 na integral é apenas

para remover termos dependentes de ∆x na n-ésima derivada de pk. No caso r = 3, a

expressão para βk fica

β0 =13

12(fi−2 − 2fi−1 + fi)

2 +1

4(fi−2 − 4fi−1 + 3fi)

2

β1 =13

12(fi−1 − 2fi + fi+1)

2 +1

4(fi−1 − fi+1)

2

β2 =13

12(fi − 2fi−1 + fi+2)

2 +1

4(3i − 4fi+1 + fi+2)

2 .

A medida da variação total tem ainda a seguinte propriedade, que provaremos na seção

4.2.2: se f é contínua no domínio de Sk e não contém pontos críticos (isto é, pontos

onde a derivada se anula), então βk = O(∆x2). Outra propriedade é que se Sk contém

20

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS2.3. LEIS DE CONSERVAÇÃO HIPERBÓLICAS

alguma descontinuidade de f então βk = O(1). É isto que nos permite utilizar a medida

da variação total como um medidor de suavidade.

2.3 Leis de conservação hiperbólicas

Os métodos que desenvolveremos no capítulo 3, como o ENO e WENO, são desenvol-

vidos para aproximar problemas da forma

∂

∂tu+

∂

∂xf(u(x, t)) = 0,

que possui apenas uma variável u. No entanto, como vimos no primeiro capítulo, muitos

dos modelos que queremos aproximar numericamente dependem de diversas variáveis,

como no caso do sistema de equações de Euler

ρ

ρv

E

t

+

ρv

ρv2 + p

v(E + p)

x

= 0.

Observamos que cada linha do sistema depende de mais de uma variável, o que impede

que os métodos que definiremos a seguir sejam utilizados de maneira trivial. O que na

verdade desejamos é que cada linha possa conter informações de apenas uma variável, isto

é, que o sistema esteja desacoplado. Felizmente os sistemas lineares hiperbólicos podem

ser diagonalizados, de modo a reescrever o mesmo sistema com novas variáveis, estas por

sua vez desacopladas. Por esta razão é que o caráter hiperbólico das leis de conservação

é tão importante para a resolução deste tipo de problema.

21

2.3. LEIS DE CONSERVAÇÃO HIPERBÓLICASCAPÍTULO 2. TÓPICOS INTRODUTÓRIOS

Definição 2. Um sistema linear da forma

ut + Aux = 0, (2.12)

com u = (u1, u2, . . . , um)T , é dito hiperbólico se a matriz A = (aij) de dimensão m ×m

é diagonalizável com autovalores reais. Se os autovalores forem todos distintos, o sitema

é dito estritamente hiperbólico.

Chamaremos u de vetor das variáveis conservadas, e denotaremos os autovalores por

λ1 ≤ λ2 ≤ . . . ≤ λm.

A matriz é diagonalizável se existem vetores não-nulos r1, r2, . . . , rm ∈ Rm, ditos os

autovetores de A, tais que

Arp = λprp, p = 1, . . . ,m

e r1, r2, . . . , rm são linearmente independentes. Neste caso, a matriz

R = [r1|r2| . . . |rm]

é não-singular, e possui uma inversa R−1. Temos então que

Λ = R−1AR e A = RΛR−1, (2.13)

22


onde

Λ =

λ1

λ2

. . .

λm

≡ diag(λ1, λ2, . . . , λm).

A importância desta representação da matriz A está no desacoplamento das variáveis.

Quando escrevemos a equação (2.12) linha por linha, isto é,

(u1)t + a11 (u1)x + . . .+ a1m (um)x = 0

(u2)t + a21 (u1)x + . . .+ a2m (um)x = 0

...

(um)t + am1 (u1)x + . . .+ amm (um)x = 0,

vemos que cada linha depende de mais de uma variável, a menos que a matriz A seja

diagonal. Exatamente por isto é que queremos fazer uma mudança de base, de modo que

matriz A seja "substituída" pela matriz diagonal Λ. Com efeito, fazendo uso das equações

(2.12) e (2.13),

ut +RΛR−1ux = 0.

Multiplicando à esquerda por R−1 e definindo w = R−1u,

0 = R−1ut +R−1RΛR−1

ux

= wt + Λwx.

Como Λ é diagonal, obtemos um sistema desacoplado de m variáveis, chamadas de variá-

23


veis características,

(wp)t + λp (wp)x = 0, p = 1, . . . ,m.

Observamos que como os autovalores λp são reais, as equações acima possuem sentido

físico.

Observação 1. Chamamos as novas variavéis w de características pois o sistema de equa-

ções acima formam um sistema de equações de advecção, cuja solução é uma combinação

linear de m ondas que viajam em velocidades características λ1, . . . , λm.

Até aqui consideramos A com coeficientes reais, mas podemos expandir toda a nossa

discussão para problemas da forma

ut + F(u)x = 0, (2.14)

com

F(u) =

f1(u)

f2(u)

...

fm(u)

,

e definindo A como o Jacobiano F′(u),

A(u) = F′(u) =

∂f1∂u1

∂f1∂u2

. . . ∂f1∂um

∂f2∂u1

∂f2∂u2

. . . ∂f2∂um

......

. . ....

∂fm∂u1

∂fm∂u2

. . . ∂fm∂um

24


Fazendo uso da regra da cadeia,

F(u)x = F′(u)ux,

obtemos

ut +A(u)ux = 0.

Como o sistema é hiperbólico, F′(u) possui m autovalores reais

λ1(u) ≤ λ2(u) ≤ . . . ≤ λm(u)

e um conjunto completo de autovetores linearmente independentes

r1(u), r2(u), . . . , rm(u).

De maneira análoga a (2.13) podemos obter

Λ(u) = R−1(u)A(u)R(u) e A(u) = R(u)Λ(u)R−1(u),

e portanto, podemos diagonalizar o sistema (2.14), obtendo novamente

wt + Λwx = 0

25


2.3.1 Equações de Euler

As equações de Euler são um sistema hiperbólico, e mostraremos aqui como podemos

realizar o desacoplamento. Relembrando mais uma vez,

ρ

ρv

E

t

+

ρv

ρv2 + p

v(E + p)

x

= 0

onde ρ é a densidade, p é a pressão, u é a velocidade e E é a energia total por unidade de

volume,

E = ρ

(

1

2u2 + e

)

,

com e sendo a energia interna específica dada pela equação de estado

e = e(ρ, p).

Para gases ideais sua expressão é dada por

e =p

(γ − 1)ρ,

e γ = cpcv

denota a razão de calor específico. Podemos também escrever a velocidade do

som como

a =

√

γp

ρ.

26


Comparando as equações de Euler com (2.12), temos

u =

ρ

ρv

E

=

u1

u2

u3

,

F(u) =

f1

f2

f3

=

ρu

ρv2 + p

v(E + p)

=

u2

12(3− γ)

u22

u1+ (γ − 1)u3

γ u2

u1u3 −

12(γ − 1)

u32

u21

,

e assim calculamos a matriz jacobiana, que é da forma

A(u) =

0 1 0

−12(γ − 3)

(

u2

u1

)2

(3− γ)(

u2

u1

)

γ − 1

−γu2u3

u21

+ (γ − 1)(

u2

u1

)3γu3

u1− 3

2(γ − 1)

(

u2

u1

)2

γ(

u2

u1

)

.

Observação 2. Observamos também por inspeção que as equações de Euler com gás ideal

satisfazem a propriedade da homogeneidade, isto é,

F(u) = A(u)u.

Usando a expressão de A(u) e o polinômio característico, resolvemos

|A−λI| = 0,

que leva a

λ1 = u− a, λ2 = u, λ3 = u+ a,

27

2.4. MÉTODOS NUMÉRICOS CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS

com os autovetores correspondentes

r1 =

1

u− a

H − ua

, r2 =

1

u

12u2

, r3 =

1

u+ a

H + ua

,

onde H é a entalpia e é definida como

H =E + p

ρ.

2.4 Métodos numéricos

Nesta seção veremos uma maneira geral de se resolver numericamente equações de leis

de conservação em uma dimensão espacial. Começaremos estudando a forma geral destas

equações e mostrando que os esquemas vistos nesta dissertação usam a mesma ferramenta

para a integração temporal, o método de Runge-Kutta TVD de ordem 3, mas utilizam

estratégias diferentes para a parte espacial.

A seguir, veremos uma maneira bem simples de se aproximar a parte espacial, chamada

de aproximação por estêncil fixo. Apesar de ele gerar oscilações na solução numérica, o

estudo deste método é importante para introduzir alguns conceitos e notações que serão

usados no capítulo 3. Depois veremos algumas propriedades que um método deve possuir

para resolver adequadamente uma lei de conservação. Por fim, veremos a causa das

oscilações nas soluções numéricas obtidas pelos esquemas de estêncil fixo.

28

CAPÍTULO 2. TÓPICOS INTRODUTÓRIOS 2.4. MÉTODOS NUMÉRICOS

2.4.1 Resolução numérica de leis de conservação

As equações de lei de conservação são da forma

∂

∂tu(x, t) +

∂

∂xf(u(x, t)) = 0. (2.15)

com condição inicial e de fronteira definidas em cada problema. Na resolução numérica

desta equação vamos definir um espaçamento ∆x e trabalhar com uma malha de pontos

x0 < x1 < . . . < xM

xi+1 = xi +∆x, i = 0, . . . ,M.

Também será útil definir os pontos de interface

xi± 12= xi ±

∆x

2, i = 0, . . . ,M.

Além disto vamos calcular a aproximação da solução u(x, t) somente em um conjunto

finito e discreto de instantes, não necessariamente uniformes,

t0 < t1 < . . . < tN .

Na discussão desta seção não nos preocuparemos com implicações da condição de fronteira

na discretização.

Fixando a equação (2.15) em xi e tn, temos

du(x, t)

dt

∣

∣

∣

∣

x=xi,tn

= −∂f(u(x, t))

∂x

∣

∣

∣

∣

x=xi,tn

29


ou, em uma notação mais simples,

du

dt

∣

∣

∣

∣

t=tn

= −∂f(u(x, tn))

∂x

∣

∣

∣

∣

x=xi

(2.16)

Os métodos vistos neste texto seguem a mesma receita para resolver (2.16). O membro

esquerdo da igualdade, a parte temporal, será resolvida usando-se o método de Runge-

Kutta de terceira ordem.

Notação 1. Chamaremos de L(u) o operador que aproxima a parte espacial de (2.16):

L(u(xi, tn)) = −∂f(u(x, tn))

∂x

∣

∣

∣

∣

x=xi

+O(∆xr) (2.17)

O método de Runge-Kutta funciona assim: a partir do tempo atual tn obtemos o

próximo estado un+1 em três passos:

u(I) = un +∆tL(un) (2.18)

u(II) =3

4un +

1

4u(I) +

1

4∆tL(u(I)) (2.19)

un+1 =1

3un +

2

3u(I) +

2

3∆tL(u(II)) (2.20)

Desta maneira, temos uma forma de obter avanços no tempo independentemente do

método que usemos para aproximar a parte espacial de (2.16). Por causa disto, no restante

desta dissertação só nos preocuparemos com a parte espacial. É nela que reside a diferença

entre um método comum e um método essencialmente não-oscilatório.

2.4.2 Aproximação por estêncil fixo

Agora veremos uma maneira simples de se aproximar o lado direito de (2.17).

30


Notação 2. Chamaremos

f(x, t) ≡ f(u(x, t)).

Quando fixarmos t e não houver confusão sobre seu valor, denotaremos simplesmente

f(x) ≡ f(x, t)

e

f(xi) ≡ fi

Assim, por (2.17), o que estamos querendo achar é um operador L(u) tal que

L(u(xi, tn)) = −f ′(xi) +O(∆xr), i = 0, . . . ,M.

Toda a informação que temos são os valores pontuais

f0, f1, . . . , fi, . . . fM .

A partir deles, é possível criar uma aproximação polinomial para f ′(xi). Para isto,

vamos usar uma função especial, chamada de fluxo numérico.

Definição 3. O fluxo numérico h(x) é a função que satisfaz

f(x) =1

∆x

∫ x+∆x2

x−∆x2

h(ξ)dξ (2.21)

Denotaremos

h(xi) ≡ hi. (2.22)

Em resumo, h é definida de tal maneira que f ′(xi) é a média de h no intervalo [xi− 12xi+ 1

2]

31


para todo i. Por causa disto, temos a interessante propriedade de que a derivada de f em

xi é exatamente (e não aproximadamente) uma diferença finita de h:

f ′(xi) =hi+ 1

2− hi− 1

2

∆x. (2.23)

Infelizmente, não sabemos nada sobre o valor hi+ 12. Mas é fácil achar aproximações para

ele. O que faremos agora é obter uma aproximação polinomial para hi+ 12, de ordem r, a

partir de r valores pontuais de f . Antes disto, iremos definir um objeto que será muito

útil de aqui por diante.

Definição 4. Um estêncil é um operador que associa um ponto a um conjunto de pontos.

Nesta dissertação, usaremos a seguinte família de estênceis:

Sr,k,i+ 12= xi+k−r+1, xi+k−r+2, . . . , xi+k

Nesta notação, r é o número de pontos do estêncil, i+ 12

é o índice do ponto de referência

e k é o índice de deslocamento lateral do estêncil. Por vezes, apenas o índice k será usado

para identificar o estêncil.

É possível encontrar os coeficientes do polinômio que aproxima hi+ 12usando os pontos

do estêncil Sr,k,i+ 12. Um pouco de cálculo usando os resultados da seção 2.1.1 (vide [5])

nos fornecem os coeficientes

cr,k,j =r∑

m=j

r∑

n=0,n 6=m

r∏

p=0,p 6=n,m

(r − k − p)

r∏

n=0,n 6=m

(m− n)(2.24)

que satisfazemr∑

j=1

cr,k,jv(xi+k−r+j) = hi+ 12+O(∆xr).

32


Definição 5. Definimos a aproximação pelo estêncil Sr,k,i+ 12, fk

r,i+ 12

, como

fkr,i+ 1

2=

r∑

j=1

cr,k,jvi+k−r+j.

O índice r indica a ordem da aproximação, k indica o deslocamento lateral e i+ 12

denota

que esta é uma aproximação para hi+ 12.

Note que como os coeficientes cr,k,j não dependem de i (por (2.24)), usamos os mesmos

coeficientes para obter fkr,i− 1

2

:

fkr,i− 1

2=

r∑

j=1

cr,k,jfi+k−r+j−1.

Agora já possuímos conhecimento o suficiente para construir o operador L(u). Primei-

ramente fixamos um deslocamento k entre 0 e r − 1, para que xi sempre esteja presente

no estêncil Sr,k,i+ 12. Depois fazemos

L(u(xi, tn) = −fkr,i+ 1

2

− fkr,i− 1

2

∆x, i = 0, . . . ,M.

Podemos verificar que L(u) assim definida satisfaz

L(u(xi, tn) = −

[

hi+ 12+O(∆xr)

]

−[

hi− 12+O(∆xr)

]

∆x(2.25)

= −f ′(xi) +O(∆xr). (2.26)

Observação 3. Pela dedução que fizemos, era para a expressão (2.25) ser apenas da

ordem de ∆xr−1 devido ao denominador ∆x. No entanto, como estamos usando o mesmo

33


deslocamento k nas duas aproximações, temos um cancelamento que faz com que

fkr,i+ 1

2− fk

r,i− 12= f ′(xi)∆x+O(∆xr+1)

Certas escolhas de k recebem nomes especiais. Quando r é par, a aproximação com

k = r2

é chamada de aproximação central. Quando r é ímpar, a aproximação com k = r−12

é chamada de aproximação upstream central.

2.4.3 O que torna um método conservativo?

Falaremos agora sobre duas peculiaridades das leis de conservação que impõem restri-

ções aos esquemas numéricos a serem usados. Não as discutiremos a fundo; para isto, [6]

é uma boa referência.

Primeiramente, quando uma equação de lei de conservação possui uma solução fraca

u, pode ser que um método numérico obtenha uma solução numérica que não é uma

aproximação de u. Para que possamos garantir que toda solução numérica obtida por

um método seja uma aproximação de uma solução analítica, devemos nos assegurar que

o método é conservativo.

Definição 6. Um esquema de diferenças finitas é dito conservativo se existe uma função

g : W −→ R Lipschitz contínua tal que:

1. O operador L(u) pode ser escrito na forma conservativa

L(u) = −g(xi+k−r+1, xi+k−r+2, . . . , xi+k)− g(xi+k−r, xi+k−r+1, . . . , xi+k−1)

∆x(2.27)

2. g é consistente com v, isto é,

g(x, x, . . . , x) ≡ f(x) (2.28)

34


Podemos verificar que os esquemas de estêncil fixo que vimos são conservativos. Basta

fazermos

g(ξ1, ξ2, . . . , ξr) =r∑

j=1

cr,k,jf(ξj).

Como os mesmos coeficientes cr,k,j são usados para obter vr,k,i+ 12

e vr,k,i− 12, é fácil ver que

L(u) = −fkr,i+ 1

2

− fkr,i− 1

2

∆x

= −g(xi+k−r+1, xi+k−r+2, . . . , xi+k)− g(xi+k−r, xi+k−r+1, . . . , xi+k−1)

∆x.

Em segundo lugar, as equações de leis de conservação não-lineares podem ter mais de

uma solução fraca para um dado problema. Entretanto, há apenas uma única solução

fraca fisicamente relevante. Trata-se da solução-limite da equação

ut + f(u)x = −κuxx

quando κ −→ 0 (isto é, quando o termo viscoso −κuxx vai a zero). Esta solução-limite

é chamada de solução fraca entropicamente correta. Uma maneira de se garantir que

um método só obtenha soluções fisicamente relevantes é estabelecer condições adicionais

ao problema, chamadas de condições de entropia, das quais apenas citamos a existência.

É possível verificar que todos os esquemas vistos nesta dissertação são entropicamente

corretos, graças à partição de fluxo de Lax-Friedrichs, conforme a indicação em [5].

2.4.4 Oscilações

O método de diferenças finitas de estêncil fixo visto aqui funciona bem quando não há

descontinuidades na solução. Caso contrário, o presente método gera oscilações numéricas.

Vejamos o seguinte exemplo, que ilustrará este fato.

35


Exemplo 1. A equação do transporte

ut + ux = 0 (2.29)

(que é uma equação da forma (2.15) com f(u) = u) com condição incial dada pela função

de salto

u0(x) =

0, se x < 0

1, se x ≥ 0

tem como solução exata a função

u(x, t) = u0(x− t).

Portanto, em t = 2 ,

u(x, 2) = u0(x− 2) =

0, se x < 2

1, se x ≥ 2(2.30)

Esta mesma equação foi resolvida numericamente usanto-se o esquema upstream central

de quinta ordem (r = 5 e k = 2). Usamos ∆x = 1128

e condição de fronteira fixa no

domínio [−1, 3]. A solução obtida em t = 2 é mostrada na figura abaixo.

(FIGURA)

Como podemos ver, nas proximidades de x = 2, onde seria a posição do choque por

(2.30), a solução apresenta oscilações significativas, com um pico de aproximadamente

0, 1 que é bem grande se comparado a ∆x ≃ 0, 008.

Estas oscilações são devidas ao fenômeno de Gibbs dos métodos espectrais. Elas

ocorrem porque, devido ao estêncil ser fixo, inevitavelmente interpolamos a função em um

domínio que contém a descontinuidade. Pelo teorema seguinte vemos que a interpolação

36


em um domínio com uma descontinuidade gera erros da ordem do tamanho do salto da

função.

Teorema 4. Seja ∆x definido como em (2.7). Seja p o maior valor tal que f ∈ Cp[x0, xk]

(considere que p = −1 se f é descontínua e que p = ∞ se f ∈ C∞[x0, xk]). Então

f(x)− p(x) =

O(

∆xp+1[f (p+1)])

, se − 1 ≤ p ≤ k − 1

O

(

∆xk+1 maxx0≤ξ≤xk

f (k+1)(ξ)

)

, se p ≥ k

para todo x ∈ [x0, xk].

A demonstração deste teorema é uma aplicação direta de séries de Taylor. Veja indi-

cação em [9].

Portanto, para que não tenhamos soluções oscilatórias devemos evitar interpolações

em estênceis descontínuos. No próximo capítulo veremos os esquemas essencialmente

não-oscilatórios e como eles conseguem evitar o uso de estênceis descontínuos.

37


38

Capítulo 3

Métodos essencialmente

não-oscilatórios

No capítulo anterior, vimos que o método para resolução numérica para leis de con-

servação unidimensional que utilizava aproximações de estêncil fixo gerava oscilações nas

vizinhanças da descontinuidade. Para lidarmos corretamente com soluções descontínuas

é necessário adotar uma outra estratégia: ao invés de usarmos sempre o mesmo estêncil

para aproximar hi+ 12

em toda a malha xi, podemos imaginar um esquema que adaptasse

o estêncil em função da localização das descontinuidades.

A adaptabilidade dos estênceis é a principal característica da classe de métodos cha-

mada de esquemas essencialmente não-oscilatórios. O nome se deve ao método deno-

minado essentially non-oscillatory schemes (ENO), proposto no artigo [1] de Harten,

Enquist, Osher e Chakravarthy. Mais tarde, outros métodos baseados no ENO surgiram.

O mais importante deles foi o esquema não-oscilatório ponderado (WENO, do original em

inglês weighted essentially non-oscillatory), desenvolvidos nos artigos [2] de Liu, Osher e

Chan e [4] de Jiang e Shu. Por sua importância e pioneirismo na história dos esque-

39

CAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS

mas essencialmente não-oscilatórios, chamaremos o WENO desenvolvido com os pesos de

Jiang-Shu [4] de WENO clássico ou WENO-JS. Recentemente, o artigo [3] de Henrick,

Aslam e Powers lançou uma nova luz sobre o WENO, com o desenvolvimento do WENO

mapeado ou simplesmente WENO-M.

Na primeira seção deste capítulo, descreveremos o método ENO. Este método escolhe

o estêncil onde a função é mais regular, dentre todos os estênceis disponíveis no domínio, e

o utiliza na aproximação de hi+ 12. Deste modo, evita-se a interpolação através de estênceis

descontínuos, o que impede que haja oscilações na solução numérica.

Na seção seguinte, veremos dois métodos, o WENO clássico e o WENO mapeado.

Mostraremos como WENO consegue aumentar a ordem de convergência em relação ao

ENO, ao fazer uma combinação das aproximações em cada um dos estênceis do domínio

ao invés de utilizar apenas um estêncil na aproximação. A combinação das aproximações

nos estênceis é feita atribuindo pesos a elas fazendo uma soma poderada. Estes pesos

são construídos de maneira a ficarem muito próximos, num domínio contínuo, dos pesos

ideais - nome dado às constantes que permitem que uma combinação de aproximações em

r estênceis consecutivos de r pontos dê uma aproximação de ordem n = 2r−1. Por outro

lado, o peso relacionado a um estêncil descontínuo deve ser bem próximo de zero, para que

se emule o comportamento do ENO em vizinhanças de descontinuidades. Deste modo,

enquanto o método ENO, que utiliza r estênceis na escolha, é de ordem r, o WENO, que

utiliza os mesmos r estênceis na aproximação, é de ordem 2r − 1.

Estudaremos o WENO clássico, e mostraremos duas falhas na construção dos pesos

que provocam perdas de ordens de convergência. A primeira delas está nas condições

impostas aos pesos, que não são suficientes para garantir a ordem ótima ao esquema em

domínios suaves. A segunda falha está no medidor de suavidade utilizado pelo WENO,

que está intimamente ligado ao valor dos pesos. Em pontos críticos, i.e., em pontos

onde a derivada se anula, os indicadores de suavidade não possuem uma propriedade

40

CAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS3.1. O MÉTODO ENO

necessária para garantir a ordem ótima naquela vizinhança. Uma análise das séries de

Taylor do valor dos medidores de suavidade nos mostra esse resultado e, além disto, indica

que, quanto maior for a ordem do ponto crítico, i.e., quanto mais derivadas se anulem

simultaneamente), piores serão os pesos obtidos.

O WENO mapeado, visto logo em seguida, introduz uma modificação ao esquema

clássico que corrige a primeira falha em todos os pontos e a segunda em boa parte dos

pontos críticos, excetuando-se aqueles em que a ordem do ponto crítico é maior ou igual

que r − 1.

Na terceira seção, veremos um novo esquema, denominado WENO-Z. Este WENO

utiliza uma fórmula diferente da fórmula clássica de Jiang-Shu. A fórmula do WENO-Z

não possui a primeira falha do esquema clássico, e garante a ordem ótima em regiões

contínuas de derivada não-nula sem a necessidade de um mapeamento, o que o torna mais

eficiente em termos computacionais. Além disto, mostraremos que o WENO-Z atribui

um peso relativamente maior aos estênceis descontínuos, acarretando um método menos

dissipativo. Os resultados preliminares para o WENO-Z se mostraram superiores aos

do WENO clássico e comparáveis ao WENO mapeado, ganhando do último com uma

vantagem estreita na maior parte dos casos.

3.1 O método ENO

No capítulo anterior vimos que um esquema de estêncil fixo de ordem r escolhe um

estêncil Sk,

Sk = xi+k−r+1, xi+k−r+2, . . . , xi+k, (3.1)

41

3.1. O MÉTODO ENOCAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS

com 0 ≤ k ≤ r − 1, e para todo xi utiliza sempre a aproximação polimonial

fki+ 1

2=

r−1∑

j=1

ckjfi−k+j, (3.2)

L(ui) = −fki+ 1

2

− fki− 1

2

∆x.

Isto implica que se acaso a função f for descontínua então certamente realizamos uma

interpolação através de um estêncil descontínuo. Como vimos, esta é a causa da oscilação

das soluções descontínuas obtidas pelos métodos de estêncil fixo.

O método ENO também se baseia nas aproximações polinomiais como em (3.2). A

diferença é que o ENO possui uma maneira de se adaptar aos estênceis em função das

descontinuidades. Partimos do princípio de que, para cada ponto xi da malha há r

estênceis possíveis que podemos escolher para obtermos L(ui) :

S0, S1, . . . , Sr−1.

Então o que fazemos é escolher o melhor dentre os r candidatos. Mas como fazemos

para saber qual é o melhor. Podemos, por exemplo, aplicar um medidor de suavidade em

cada um dos estênceis e escolher o que proporcionar o menor valor absoluto no medidor.

Isto é quase o que o método ENO faz.

Na verdade, o ENO faz o seguinte: começamos com um estêncil de um ponto, S1 =

(xi). Daí usamos o medidor de suavidade nos estênceis (xi−1, xi) e (xi, xi+1). Se

|M(xi−1, xi)| ≤ |M(xi, xi+1)| ,

42

CAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS3.1. O MÉTODO ENO

então definimos o estêncil de 2 pontos como S2 = (xi=1, xi). Caso contrário, se

|M(xi−1, xi)| > |M(xi, xi+1)| ,

fazemos S2 = (xi, xi+1). Repetimos o processo até o j-ésimo passo, quando temos um

estêncil de j pontos Sj = (xi−l, . . . , xi+m). Daí obtemos um estêncil d n+1 pontos fazendo

|M(xi−l−1, . . . , xi+m)| ≤ |M(xi−l, . . . , xi+m+1)| =⇒ Sj+1 = (xi−l−1, . . . , xi+m)

|M(xi−l−1, . . . , xi+m)| > |M(xi−l, . . . , xi+m+1)| =⇒ Sj+1 = (xi−l, . . . , xi+m+1)

(3.3)

Prosseguiremos assim até que tenhamos construído o estêncil Sr. Depois usamos a apro-

ximação do estêncil Sr para obtermos fi+ 12.

O ENO original utiliza como medidor de suavidade as diferenças divididas por causa

de sua eficiência computacional (não entraremos em detalhes sobre este medidor), mas

em tese qualquer outro medidor poderia ser utilizado.

Vamos agora mostrar algumas propriedades do método ENO. Considere Ii = [xi− 12, xi+ 1

2]

e o polinômio definido como

pi(x) =r∑

j=1

cki,jf(x+ (j + ki − r)∆x),

onde ki é tal que pi(xi) = fi+ 12. Se f é contínua no intervalo Ii, podemos dizer que

pi(x) = h(x) +O(∆xr), x ∈ Ii.

Caso contrário, se f for descontínua no intervalo Ii, podemos dizer que pi(x) é monótona

43

3.1. O MÉTODO ENOCAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS

em Ii. Além, disso, pi(x) é de variação total limitada, i.e.,

TV [pi] = supP

∑

j|xj ,xj+1∈P

|pi(xj+1)− pi(xj)| < ∞

onde P = xi− 12= x0 < x1 < . . . < xN = xi+ 1

2 é uma partição de Ii. Em outras palavras,

existe uma função z(x) que satisfaz

z(x) = pi(x) +O(∆xr), x ∈ Ii

em toda célula Ii tal que

TV [z] ≤ TV [f ].

Em outras palavras, a propriedade acima nos diz que o esquema ENO é de variação total

limitada, o que justifica a ausência de oscilações significativas nas suas soluções e é a razão

do nome "essencialmente não-oscilatório".

Procuramos investigar o caráter conservativo do ENO, mas uma dúvida persiste. Por

construção, o ENO pode utilizar estênceis diferentes para obter fi+ 12

e fi− 12. Neste caso,

o esquema ainda é conservativo? Se for, isto quer dizer que para qualquer 0 ≤ k ≤ r − 1

existe uma função g tal que

fi+ 12− fi− 1

2

∆x=

g (xi−r+1, . . . , xi+r−1)− g (xi−r, . . . , xi+r−2)

∆x

e

g(x, x, . . . , x) = v(x).

Embora tenhamos motivos para acreditar que o ENO seja conservativo (até porque

não conseguimos até agora construir nenhum contra-exemplo), não encontramos nenhuma

44

CAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS3.2. ESQUEMAS WENO

referência contendo a prova desta afirmação.

3.2 Esquemas WENO

O método ENO tem como característica a boa adaptação em relação às descontinui-

dades, mas ainda possui alguns pontos que precisam ser melhorados. O principal ponto

em questão é que na escolha do melhor estêncil são levados em conta r candidatos

S0, S1, . . . , Sr−1,

mas apenas um estêncil é escolhido. Em termos de pontos, temos um conjunto total de

2r − 1 pontos no domínio

xi+k−r+1, xi+k−r+2, . . . , xi+k,

mas usamos apenas r pontos na aproximação. Claramente existe um desperdício de

informação, já que se fossem utilizados todos os 2r − 1 pontos poderíamos obter uma

aproximação de ordem 2r − 1 ao invés da ordem r que obtemos com o ENO.

Algumas modificações foram desenvolvidas para realizar melhorias no ENO, mas den-

tre todas as opções o método WENO foi a mais bem-sucedida. Atualmente há várias

versões do WENO, mas todas elas seguem a mesma ideia básica: em vez de usarmos ape-

nas um estêncil Sk para obtermos a aproximação fi+ 12, vamos fazer uma combinação das

aproximações fki+ 1

2

, atribuindo um peso ωk a cada uma delas:

fi+ 12=

r−1∑

k=0

ωkfki+ 1

2(3.4)

45

3.2. ESQUEMAS WENOCAPÍTULO 3. MÉTODOS ESSENCIALMENTE NÃO-OSCILATÓRIOS

Os pesos devem satisfazer a algumas condições básicas. Primeiramente, por consis-

tência e estabilidade,

Condição 1 (convexidade). Os pesos devem satisfazer às condições de convexidade

ωk ≥ 0,r−1∑

k=0

ωk = 1.

Além disto, para que este método seja tão adaptável a descontinuidades quanto o ENO,

gostaríamos de atribuir peso praticamente nulo às aproximações obtidas em estênceis com

descontinuidade. Isto evitaria que usássemos na soma (3.4) uma parcela que tivesse sido

obtida sob a influência do fenômeno de Gibbs. Portanto,

Condição 2 (emulação do ENO). Se o k-ésimo estêncil é descontínuo (e há pelo menos

um estêncil contínuo), então

ωk ≈ 0.

A última condição é a que garante um ganho em termos de ordem de convergência.

Estipulamos que

Condição 3 (optimalidade). Em regiões contínuas,

fi+ 12− fi− 1

2

∆x=

r−1∑

k=0

ω+k f

ki+ 1

2

−r−1∑

k=0

ω−k f

ki− 1

2

∆x= f ′(xi) +O(∆x2r−1),

onde ω±k denota os pesos usados na construção de fi± 1

2.

A construção dos pesos ωk geralmente passa pelo uso dos chamados pesos ideais, nome

dado às constantes que satisfazem às condições (1) e (3). Por causa disto, os pesos ideais

serão o assunto da seção a seguir.

46


3.2.1 Pesos ideais e comparação com upstream central

Para satisfazer a condição (3), precisamos utilizar informações de todos os 2r − 1

pontos. Em outras palavras, queremos que em regiões contínuas, a escolha dos pesos

ideais seja feita de forma que ela exerça o mesmo papel que um interpolador de grau

2r − 1, em um estêncil maior Sn = xi+k−r+1, xi+k−r+2, . . . , xi+k, formado pelos 2r − 1

pontos em questão. Chamaremos este interpolador de Fi+ 12.

Afirmamos que sempre é possível realizar uma combinação das aproximações f 0i+ 1

2

, . . . , f r−1i+ 1

2

que resulta na aproximação upstream central Fi+ 12. Em outras palavras, existem constan-

tes dk, chamadas de pesos ideais, que satisfazem

r−1∑

k=0

dkfki+ 1

2= Fi+ 1

2. (3.5)

Observação 4. Não encontramos nenhuma prova deste fato nas referências, embora ele

seja largamente utilizado. De qualquer forma, uma demonstração não parece difícil de se

obter.

Observação 5. Os pesos ideais dk dependem da ordem do método, ou seja, dk = dr,k. No

entanto, por praticidade, omitiremos a variável r na maioria dos casos.

Como dito anteriormente, supondo as condições anteriores como verdadeiras, temos

que os pesos ideais serviriam perfeitamente para o WENO numa região contínua. Em

outras palavras, deseja-se que os pesos do WENO estejam bem próximos dos pesos ideais

numa região contínua para que a aproximação do WENO fique próxima da aproximação

upstream central. Veremos algumas maneiras de se garantir esta proximidade na próxima

seção e seguintes.

Vamos mostrar como se calcula os pesos ideais de ordem r = 3. Os pesos de outras

ordens apenas explicitaremos. Para tal, considere o seguinte problema:

47


Problema 2. Encontre os pesos ideais de ordem 3, isto é, as constantes d0, d1 e d2 tais

que

d0f0i+ 1

2+ d1f

1i+ 1

2+ d2f

2i+ 1

2= Fi+ 1

2

Usando a definição de aproximação

fki+ 1

2=

r∑

j=1

ckjfi+k−r+j,

obtemos

f 0i+ 1

2=

1

6(2fi−2 − 7fi−1 + 11fi) ,

f 1i+ 1

2=

1

6(−fi−1 + 5fi + 2fi+1) ,

f 2i+ 1

2=

1

6(2fi + 5fi+1 − fi+2) ,

Fi+ 12

=1

60(2fi−2 − 13fi−1 + 47fi + 27fi+1 − 3fi+2) .

Deste modo, temos o seguinte sistema:

1

6

2 0 0

−7 −1 0

11 5 2

0 2 5

0 0 −1

d0

d1

d2

=1

60

2

−13

47

27

−3

Apesar de este ser um sistema sobredeterminado, ele possui uma solução, que é

d0 =1

10, d1 =

6

10, d2 =

3

10.

48


A tabela a seguir mostra o valor dos pesos ideais para r = 1, 2, 3, 4, 5 e 6 .

Pesos ideais para r = 1, 2, 3, 4, 5 e 6

r dr,0 dr,1 dr,2 dr,3 dr,4 dr,5

1 1

2 13

23

3 110

610

310

4 135

1235

1835

435

5 1126

20126

60126

40126

5126

6 1462

30462

150462

200462

75462

6462

Agora estamos prontos para estudar os esquemas WENO.

3.2.2 WENO clássico: os pesos de Jiang-Shu

Na versão original do WENO [2], a condição de optimalidade ainda não havia estipu-

lada daquela maneira. Em vez dela, tínhamos

fi+ 12=

r−1∑

k=0

ωkfki+ 1

2= hi+ 1

2+O(∆xr+1).

Note o termo O(∆xr+1) nesta expressão. Isto quer dizer que, no máximo, o WENO

original obtinha apenas ordem máxima de r + 1, um ganho de somente uma ordem em

relação ao ENO.

Os primeiros a notar que o WENO poderia obter um ganho maior, de ordem até 2r−1,

foram G.-S. Jiang e C.-W. Shu. No artigo [4], eles propõem uma outra maneira de se cal-

cular os pesos, baseados num novo medidor de suavidade: a medida de variação total. Os

pesos são obtidos como uma função racional dos medidores de suavidade. O uso de uma

49


fórmula racional para obter os pesos - e não um algoritmo que usasse estruturas condici-

onais, por exemplo - deveu-se a questões de eficiência computacional.[5, p.16, §3 e p.18,

§4]. Na época (1996) as estruturas condicionais eram ineficientes em supercomputadores

com o CRAY, e, por isto, o uso de tais estruturas era evitado a qualquer custo.

Vale também ressaltar que, ao invés de usar a condição de optmalidade (3), foi utilizada

outra, mais fraca:

Condição 4 (Jiang-Shu). Em regiões contínuas,

fi+ 12=

r−1∑

k=0

ωkfki+ 1

2= hi+ 1

2+O(∆x2r−1).

Note que esta condição não é suficiente para garantir a optimalidade do esquema:

supondo esta condição como válida, em geral temos apenas

fi+ 12− fi− 1

2

∆x=

[

hi+ 12+O(∆x2r−1)

]

−[

hi− 12+O(∆x2r−1)

]

∆x= f ′(xi) +O(∆x2r−2).

Mesmo assim, os pesos de Jiang-Shu merecem um estudo mais profundo. Primeiro, porque

a condição de optimalidade é satisfeita às vezes, especialmente com ∆x pequeno. Segundo,

porque os pesos de Jiang-Shu são a base da maioria dos esquemas WENO que vieram

depois.

Vejamos então os pesos de Jiang-Shu. A fórmula dos pesos é

ωk =αk

∑r−1j=0 αj

, αk =dk

(βk + ǫ)p(3.6)

onde ǫ > 0 é um valor teoricamente utilizado apenas para evitar que o denominador se

anule na segunda equação, p é uma potência inteira e βk é a medida de variação total. Os

valores mais geralmente utilizados para estes parâmetros são p = 2 e ǫ = 10−6.

50


Notação 3. Onde houver confusão sobre o tipo do peso, denotaremos os pesos de Jiang-

Shu como ω(JS)k .

Vejamos como os pesos de Jiang-Shu satisfazem às condições do WENO. Trivialmente

a condição de convexidade é satisfeita, já que

r−1∑

k=0

ωk =r−1∑

k=0

αk∑r−1

k=0 αk

=

∑r−1k=0 αk

∑r−1k=0 αk

= 1

e

dk ≥ 0, (βk + ǫ)p > 0.

A condição (2) também é facilmente verificável. Suponha que o estêncil m é descontí-

nuo e o estêncil n é contínuo. Afirmamos que βm = O(1) e βn = O(∆x2). Não provaremos

a primeira afirmação, enquanto que a segunda será provada na seção 4.2.2. Então

ωm =

dm(βm+ǫ)p

dm(βm+ǫ)p

+ dn(βn+ǫ)p

+∑r−1

k=0k 6=m,n

dk(βk+ǫ)p

=dm

dm + O(1)(O(∆x2)+ǫ)p

+∑r−1

k=0k 6=m,n

O(1)(βk+ǫ)p

≤dm

dm +O( 1Cp ) +O(1)

= O(Cp), C = max∆x2, ǫ

A condição de Jiang-Shu (4) pode ser satisfeita da seguinte maneira: usando-se a

51


propriedade dos pesos ideais (3.5) e a condição (1), temos

r−1∑

k=0

ωkfki+ 1

2=

r−1∑

k=0

dkfki+ 1

2+

r−1∑

k=0

(ω±k − dk)f

ki+ 1

2

=r−1∑

k=0

(ωk − dk)[

hi+ 12+O(∆xr)

]

+[

hi+ 12+O(∆x2r−1)

]

=

[

hi+ 12

r−1∑

k=0

(ωk − dk) +r−1∑

k=0

(ωk − dk)O(∆xr)

]

+[

hi+ 12+O(∆x2r−1)

]

=r−1∑

k=0

(ωk − dk)O(∆xr) + hi+ 12+O(∆x2r−1) (3.7)

Assim, se tivermos

ωk − dk = O(∆xr−1), 0 ≤ k ≤ r − 1,

a condição (4) é imediatamente satisfeita.

O próximo teorema nos mostrará que o medidor de suavidade baseado na medida de

variação total faz com que os pesos de Jiang-Shu atendam à condição acima:

Teorema 5. Seja β0, . . . , βr−1 os medidores de suavidade baseados na medida de variação

total e d0, . . . , dr−1os pesos ideais de ordem r. Considere também os pesos de Jiang-Shu

ω0, . . . , ωr−1 dados pela fórmula

ωk =

dk(βk+ǫ)p

∑r−1k=0

dj(βj+ǫ)p

, k = 0, . . . , r − 1

com 0 < ǫ ≪ 1 e p inteiro maior ou igual a 1. Agora suopnha que para qualquer 0 ≤ k ≤

r − 1, possamos escrever βk na forma

βk = D(1 +O(∆xq))

52


onde D é uma constante não-nula independente de j e q ≥ 1. Então

ωk − dk = O(∆xq).

Demonstração. Usando a equação (3.6) e tomando ǫ = 0 para efeito de simplificação,

obtemos

αk =dk

(D(1 +O(∆xq)))p=

dkDp

(

1 +O(∆xq) +O(∆x2q) + ...)p

=dkDp

(1 +O(∆xq)).

Somando todos os termos,r−1∑

k=0

αk =1

Dp(1 +O(∆xq)),

levando em conta que os pesos ideais somam 1, pela condição (1). Com isso,

ωk =dkDp (1 +O(∆xq))1Dp (1 +O(∆xq))

= dk(1 +O(∆xq)).

Consequentemente,

ωk − dk = O(∆xq),

como queríamos.

Corolário 1. Se βk = D(1 +O(∆xr−1)), a condição (4) é imediatamente satisfeita.

Neste momento, vamos deter a nossa análise ao medidor de suavidade βk para verifi-

carmos a condição (4). Consideremos a análise de βk no caso r = 3. A expressão para βk

53


é

β0 =13

12(fi−2 − 2fi−1 + fi)

2 +1

4(fi−2 − 4fi−1 + 3fi)

2

β1 =13

12(fi−1 − 2fi + fi+1)

2 +1

4(fi−1 − fi+1)

2

β2 =13

12(fi − 2fi−1 + fi+2)

2 +1

4(3i − 4fi+1 + fi+2)

2 .

Expandindo as equações acima em séries de Taylor em torno de xi, temos

β0 = (f ′∆x)2 +

(

13

12f ′′2 −

2

3f ′f ′′′

)

∆x4 +

(

−13

6f ′′f ′′′ +

1

2f ′f (4)

)

∆x5 +O(∆x6),(3.8)

β1 = (f ′∆x)2 +

(

13

12f ′′2 +

1

3f ′f ′′′

)

∆x4 +O(∆x6),

β2 = (f ′∆x)2 +

(

13

12f ′′2 −

2

3f ′f ′′′

)

∆x4 +

(

13

6f ′′f ′′′ −

1

2f ′f (4)

)

∆x5 +O(∆x6).

Vamos supor inicialmente que f ′ 6= 0. Se fizermos D = (f ′∆x)2, podemos ver, usando as

equações acima, que

β0 = D

(

1 +

(

1312f ′′2 − 2

3f ′f ′′′

)

∆x4 +O(∆x5)

(f ′∆x)2

)

= D(1 +O(∆x2),

β1 = D

(

1 +

(

1312f ′′2 + 1

3f ′f ′′′

)

∆x4 +O(∆x6)

(f ′∆x)2

)

= D(1 +O(∆x2),

β2 = D

(

1 +

(

1312f ′′2 − 2

3f ′f ′′′

)

∆x4 +O(∆x5)

(f ′∆x)2

)

= D(1 +O(∆x2).

Portanto, temos βk = D(1 +O(∆xr−1)) e a condição (4) é imediatamente satisfeita.

Agora vamos supor que f ′ = 0 e f ′′ 6= 0. Expandindo novamente em séries de Taylor,

temos agora

54


β0 =13

12f ′′2∆x4 −

13

6f ′′f ′′′∆x5 +O(∆x6),

β1 =13

12f ′′2∆x4 +O(∆x6),

β2 =13

12f ′′2∆x4 +

13

6f ′′f ′′′∆x5 +O(∆x6).

Se fizermos D assumir o valor do único termo que não varia nestas três expansões, ou

seja,

D =13

12f ′′2∆x4

o melhor que conseguimos é

β0 = D

(

1−136f ′′f ′′′∆x5 +O(∆x6)

1312f ′′2∆x4

)

= D(1 +O(∆x)),

β1 = D

(

1 +O(∆x6)1312f ′′2∆x4

)

= D(1 +O(∆x2)) = D(1 +O(∆x)),

β2 = D

(

1 +136f ′′f ′′′∆x5 +O(∆x6)

1312f ′′2∆x4

)

= D(1 +O(∆x)),

a nçao ser que f ′′′ = 0 o que nos daria

βk = D(1 +O(∆x2)), k = 0, 1, 2

Deste modo, a condição (4) só é assegurada quando f ′′′ se anula. Do contrário, pelo

teorema 5 temos apenas

ωk − dk = O(∆x), k = 0, 1, 2

55


e assim a condição (4) não é necessariamente satisfeita. Na verdade, por (3.7), teríamos,

r−1∑

k=0

ωkfki+ 1

2= hi+ 1

2+O(∆x2r−2),

o que já representa uma queda na ordem de convergência do WENO.

Jiang e Shu não perceberam as duas falhas que o WENO com seus pesos possui, que

causam queda na ordem ótima 2r − 1: a condição (4), que não garante a optimalidade;

e os seus pesos que, em geral, não satisfazem à própria condição (4) em pontos críticos.

As falhas não foram percebidas porque o valor utilizado para ǫ (10−6) é relativamente

alto em comparação ao valor de βk em regiões contínuas, fazendo com que, por (3.6), os

pesos ficassem próximos dos pesos ideais o suficiente para que ambas as falhas não fossem

notadas. Além disto, a série de Taylor de βk não havia sido feita até os termos mais

importantes (faltaram os termos da ordem de ∆x5) e a função utilizada para teste era

uma senóide com descontinuidade, na qual f ′ e f ′′′ se anulavam simultaneamente [4, pp.

9-10]. Por isto não foi verificado que os pesos não satisfaziam à condição de Jiang-Shu

em pontos críticos mais gerais.

As duas falhas de projeto só foram documentadas na literatura quase uma década

depois do trabalho de Jiang e Shu, no artigo de A. Henrick, T. Aslam & J. Powers [3].

Este artigo estabelece uma relação empírica entre a ordem do WENO na vizinhança de

um ponto crítico e a chamada ordem do ponto crítico, que indica quantas derivadas se

anulam simultaneamente. O estudo desta relação será o assunto da próxima seção.

3.2.3 Pontos críticos

Em [3] é sugerida uma relação direta entre o comportamento de βk nas vizinhanças

de pontos críticos com o número de derivadas que se anulam simultaneamente. Para ser

56


mais preciso, vamos usar a seguinte definição:

Definição 7. A ordem ncp de um ponto crítico xi é o valor tal que

f ′(xi) = . . . = f (ncp)(xi) = 0, f (ncp+1)(xi) 6= 0. (3.9)

Caso todas as derivadas se anulem em xi, dizemos que este é um ponto crítico de ordem

infinita.

Através de uma série de experimentos numéricos, foi estabelecida uma relação entre

a ordem do ponto crítico e a ordem do WENO clássico na vizinhança do ponto crítico.

Apesar de não ter sido fornecida uma prova para o caso geral, séries de Taylor para βk

de ordem 3 até 6 confirmam o resultado empírico obtido em [3]. Algumas das expressões

destas séries podem ser encontradas na seção 4.3.

Nestes termos, vamos estabelecer uma conjectura prevista, de um modo um pouco

diferente, por Henrick, Aslam e Powers em [3].

Conjectura 1 (Henrick-Aslam-Powers). Seja f uma função contínua por partes com

um ponto crítico xi de ordem ncp e βk o indicador de suavidade de ordem r de f , com

espaçamento ∆x. Então βk avaliada em uma vizinhança de xi assume a forma

βk = D(1 +O(∆xq)), k = 0, . . . , r − 1,

com D independente de k e com a potência q dada por

q = max(r − 1− ncp, 0).

57


Recorrendo a (3.7)

fi+ 12=

r−1∑

k=0

ωkfki+ 1

2= hi+ 1

2+

r−1∑

k=0

(ωk − dk) fki+ 1

2+O(∆x2r−1)

e ao teorema 5, temos o seguinte corolário da conjectura 1 :

Corolário 2. Seja xi um ponto crítico de f com ordem ncp. Considere

fj+ 12=

r−1∑

k=0

ω(JS)k fk

j+ 12.

Então, para um xj vizinho a xi, a ordem rc de convergência do WENO clássico

fj+ 12− fj− 1

2

∆x= f ′(xj) +O(∆xrc)

é dado por

rc = max(2r − 2− nc, nc).

Demonstração. Supomos a conjectura 1 como sendo válida. Usando o teorema 5, temos

ω(JS)k = dk +O(∆xp), p = max(r − 1− nc, 0).

Vamos dividir o problema em dois casos:

Caso 1 (Caso 1). nc < r − 1

Neste caso, p = r − 1− nc e 2r − 2− nc > nc. Por (3.7),

fj+ 12=

r−1∑

k=0

ω(JS)k fk

j+ 12= hj+ 1

2+O(∆xr+p).

58


Logo,

fj+ 12− fj− 1

2

∆x=

[

hj+ 12+O(∆xr+p)

]

−[

hj+ 12+O(∆xr+p)

]

∆x

= f ′(xj) +O(∆xr+p−1)

= f ′(xj) +O(∆x2r−2−nc).

Caso 2 (Caso 2). nc ≥ r − 1

Agora temos p = 0 e nc ≥ 2r − 2 − nc. Como as nc primeiras derivadas de v se anulam

simultaneamente, a aproximação polinomial fkj+ 1

2

é de ordem nc + 1 ao invés de simples-

mente r:

fkj+ 1

2= hj+ 1

2+O(∆xnc+1), k = 0, . . . , r.

Com isto, não importando o valor dos pesos, o resultado da combinação é

fj+ 12

=r−1∑

k=0

ω(JS)k fk

j+ 12

=r−1∑

k=0

ω(JS)k

[

hj+ 12+O(∆xnc+1)

]

=[

hj+ 12+O(∆xnc+1)

]

r−1∑

k=0

ω(JS)k

=[

hj+ 12+O(∆xnc+1)

]

.

Assim,fj+ 1

2− fj− 1

2

∆x= f ′(xj) +O(∆xnc),

como queríamos.

59


O corolário nos informa que há uma ordem do ponto crítico que faz o WENO clássico

gerar os piores resultados. Isto acontece quando

nc = r − 1,

o que dá

rc = r − 1.

Na seção seguinte, veremos uma maneira de se obter resultados melhores em pontos

críticos.

3.2.4 WENO mapeado: corrigindo algumas falhas

Além do estudo sobre as falhas do WENO clássico, [3] contém também uma correção

para elas, sob a forma de uma função mapeadora. O mapeamento nos pesos de Jiang-Shu,

sob certas condições, os aproximam dos pesos ideais, de modo a satisfazer a condição de

optimalidade 3

Em [3] primeiramente é sugerido um critério suficiente aos pesos de modo a satisfazer

a condição 3. Expandindo a expressão

fi+ 12− fi− 1

2

∆x(3.10)

60


temos

fi+ 12− fi− 1

2

∆x=

r−1∑

k=0

ω+k f

ki+ 1

2

−r−1∑

k=0

ω−k f

ki− 1

2

∆x

=

r−1∑

k=0

d+k fki+ 1

2

−r−1∑

k=0

d−k fki− 1

2

∆x+

r−1∑

k=0

(ω+k − dk)f

ki+ 1

2

−r−1∑

k=0

(ω−k − dk)f

ki− 1

2

∆x

=Fi+ 1

2− Fi− 1

2

∆x+

r−1∑

k=0

(ω+k − dk)

[

hi+ 12+O(∆xr)

]

−r−1∑

k=0

(ω−k − dk)

[

hi+ 12+O(∆xr)

]

∆x

= f ′(xi) +O(∆x2r−1) +

r−1∑

k=0

(ω+k − dk)O(∆xr)−

r−1∑

k=0

(ω−k − dk)O(∆xr)

∆x

Portanto, uma condição suficiente para garantir a optimalidade é que

ω±k = dk +O(∆xr), k = 0, . . . , r − 1. (3.11)

A partir daí foi construída uma função mapeadora que, aplicada aos pesos de Jiang-Shu,

os faz satisfazer (3.11).

A função mapeadora gk(ω) é definida como

gk(ω) =ω (dk − d2k − 3dkω + ω2)

d2k + ω (1− 2dk), k = 0, . . . , r − 1.

Esta função tem as seguintes propriedades:

· gk é monótona crescente.

· 0 ≤ gk(ω) ≤ 1, gk(0) = 0 e gk(1) = 1.

· gk(ω) ≈ 0 se ω ≈ 0, gk(ω) ≈ 0 se ω ≈ 0.

61


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

k=0k=1k=2Identidade

Figura 3.1: Funções de mapeamento usados no WENO mapeado de ordem r = 3 parak = 0, 1 e 2

· gk(dk) = dk, g′k(dk) = g′′k(dk) = 0.

· gk(ωk) = dk +O(∆x3p) se ωk = dk +O(∆xp).

Vamos provar a última propriedade. Suponha que ωk = dk +O(∆xp). Então

gk(ωk) = gk(dk) + g′k(dk) (ωk − dk) +g′′k(dk)

2(ωk − dk)

2 +g′′′k (dk)

6(ωk − dk)

3 + . . .

= dk +g′′′k (dk)

6(ωk − dk)

3 + . . .

= dk +O(∆x3).

Vejamos como podemos usar a função gk para melhorar os pesos de Jiang-Shu. Suponha

que usamos (3.6) para obtermos os pesos ω(JS)0 , . . . , ω

(JS)r−1 . Se o estêncil Sk é descontínuo

então, como vimos, ω(JS)k ≈ 0. Assim, gk(ω

(JS)k ) ≈ 0. Por outro lado, se todos os estênceis

forem contínuos e se

ω(JS)k = dk +O(∆xp)

62


para algum p ≥ 1, k = 0, . . . , r− 1, então gk aproxima ainda mais o peso ω(JS)k em direção

a dk, fazendo

gk(ω(JS)k ) = dk +O(∆x3p),

para k = 0, . . . , r − 1.

A condição (3.11) pode ser satisfeita através de repetidas aplicações de Jiang-Shu. Se

p ≥ 1. então temos

ω(JS)k = dk +O(∆x)

ω1k = gk(ω

(JS)k ) = dk +O(∆x3)

ω2k = gk(ω

1k) = dk +O(∆x9)

...

ω⌈log3 r⌉k = gk(ω

⌈log3 r⌉k ) = dk +O(∆xr).

Portanto, ⌈log3 r⌉ aplicações sucessivas do mapeamento gk são capazes de corrigir a ordem

do esquema. Como geralmente as implementações usam r ≤ 6, duas aplicações de gk já

bastam para corrigir o problema, enquanto que para r = 3 apenas uma aplicação já é

suficiente.

Quando p = 0, contudo, o mapeamento não é capaz de fazer nenhuma correção, pois

gk(dk +O(1)) = dk +O(1).

Vimos na seção 3.2.3 que p = 0 na vizinhança de qualquer ponto crítico de ordem nc ≥

r − 1. Portanto, o mapeamento nos pesos não é capaz de garantir a optimalidade nestas

regiões contínuas. Isto não chega a ser um grande problema, já que não é comum em

aplicações práticas lidar com funções com pontos críticos de ordens altas.

63


O WENO mapeado apresenta melhores resultados que o WENO clássico mesmo quando

não há a influência de pontos críticos. No próximo capítulo veremos um novo método pro-

posto, chamado WENO-Z.

64

Capítulo 4

O método WENO-Z

No capítulo anterior apresentamos os esquemas essencialmente não oscilatórios (ENO)

e os ponderados (WENO). Vimos que enquanto o método ENO pode alcançar ordem de

convergência de até r, o método WENO, por utilizar informações de todos os r estênceis

candidatos, pode chegar a ordem 2r − 1 de convergência. Mas, para isso, tornou-se

necessário o uso de um mapeamento, pois o WENO clássico proposto por [4] não conseguia

alcançar a optimalidade. Este mapeamento foi proposto em [3] (dando o nome de WENO

mapeado), e em casos onde a ordem do ponto crítico é menor do que r−1 (o que abrange

uma grande quantidade dos problemas que resolvemos), é sempre possível alcançar a

optimalidade. No entanto, este mapeamento é computacionalmente caro, e se diversas

iterações forem realizadas, o esquema pode se tornar muito lento.

Neste capítulo, introduziremos uma nova maneira de se obter os pesos do WENO, sem

a necessidade de um mapeamento, e que também garante a optimalidade em determinadas

condições. A ideia fundamental é que utilizamos um medidor de suavidade global, que

utiliza informações de todos os 2r − 1 pontos também no medidor de suavidade. Este

medidor de suavidade global, chamado de τn, nada mais é do que uma combinação linear

dos indicadores de suavidade βk, construídos usando a medida de variação total.

65

4.1. UMA NOVA FÓRMULA PARA OS PESOS DO WENOCAPÍTULO 4. O MÉTODO WENO-Z

Na primeira seção, mostraremos o WENO-Z para ordem r = 3 e os resultados en-

contrados, e sua generalização para ordem r. Na segunda seção mostraremos alguns

resultados auxiliares sobre a aproximação polinomial de hi+ 12, denotada por fi+ 1

2, e os

indicadores de suavidade βk, Na terceira seção mostraremos uma fórmula geral para o

medidor de suavidade global τn que garante a optimalidade em situações em que não há

pontos críticos. Na quarta seção mostraremos, dentre as possíveis escolhas para se cons-

truir o medidor τn, qual é a que possui ordem máxima, ou seja, aquela que perde menos

ordem de convergência em problemas com pontos críticos. Veremos que em determinados

casos (com r relativamente grande), mesmo que haja a presença de pontos críticos, a opt-

malidade pode ser alcançada. Por fim, na última seção veremos uma fórmula alternativa

dos medidores de suavidade global que permite alcançar a optimalidade em casos em que

antes não era possível sem o mapeamento proposto no WENO mapeado.

4.1 Uma nova fórmula para os pesos do WENO

Nesta seção introduziremos um novo conjunto de pesos WENO. Para facilitar o en-

tendimento, iremos considerar inicialmente o caso de ordem r = 3, ou n = 2r − 1 = 5.

Definição 8. Vamos definir o medidor de suavidade global τ5 como sendo

τ5 = |β0 − β2| .

Pelas expansões em séries de Taylor de β0 e β2,

β0 = (f ′∆x)2 +

(

13

12f ′′2 −

2

3f ′f ′′′

)

∆x4 +

(

−13

6f ′′f ′′′ +

1

2f ′f (4)

)

∆x5 +O(∆x6),

β2 = (f ′∆x)2 +

(

13

12f ′′2 −

2

3f ′f ′′′

)

∆x4 +

(

13

6f ′′f ′′′ −

1

2f ′f (4)

)

∆x5 +O(∆x6),

66

CAPÍTULO 4. O MÉTODO WENO-Z4.1. UMA NOVA FÓRMULA PARA OS PESOS DO WENO

podemos deduzir alguns resultados: se os estênceis S0, S1 e S2 são contínuos e afastados

de qualquer ponto crítico, então

τ5 = O(∆x5).

Por outro lado, se um dentre os estênceis acima descritos for descontínuo, então

τ5 = O(1).

Propomos uma nova fórmula para os pesos que se utiliza da informação global dos

estênceis com o objetivo de obter uma melhor aproximação. A fórmula é a seguinte:

ωzk =

αzk

∑2j=0 α

zj

, αzk =

dkβzk

= dk

(

1 +τ5

βk + ǫ

)

, k = 0, 1, 2. (4.1)

Esta nova fórmula pode ser vista como uma modificação da fórmula de Jiang-Shu, com

p = 1 e

βzk =

βk

βk + τ5 + ǫ.

Ao contrário do WENO clássico, em que é crucial uma escolha relativamente alta para

ǫ (o valor sugerido é de ǫ = 10−6), no WENO-Z podemos escolher um valor baixo, como

por exemplo, ǫ = 10−40.

Uma vantagem desta fórmula é que os pesos satisfazem imediatamente à condição

de optimalidade 3 (longe de pontos críticos), sem a necessidade de um mapeamento.

Recordando o que vimos no capítulo anterior (mais precisamente na equação (3.11)) uma

condição sobre os pesos suficiente para garantir a optimalidade é que

ωk = dk +O(∆x3).

67


Isto pode ser verificado diretamente pela fórmula (usando por praticidade ǫ = 0)

ωzk =

dk

(

1 + τ5βk

)

∑r−1j=0 dj

(

1 + τ5βj

)

=dk

(

1 + O(∆x5)O(∆x2)

)

∑r−1j=0 dj

(

1 + O(∆x5)O(∆x2)

)

=dk (1 +O(∆x3))

(1 +O(∆x3))∑r−1

j=0 dj

=dk (1 +O(∆x3))

(1 +O(∆x3))

= dk(

1 +O(∆x3)) (

1 +O(∆x3))

= dk +O(∆x3).

Note que estamos utilizando que βk = O(∆x2), fato que ainda não provamos. Podemos

também verificar que a emulação do ENO 2 é preservada. Supondo-se o estêncil Sm

descontíuo e o estêncil Sn contínuo, temos

ωzm =

dm

(

1 + τ5βk

)

dm

(

1 + τ5βm

)

+ dn

(

1 + τ5βn

)

+r−1∑

j=0j 6=m,n

dj

(

1 + τ5βj

)

=dm

(

1 + O(1)O(1)

)

dm

(

1 + O(1)O(1)

)

+ dn

(

1 + O(1)O(∆x2)

)

+r−1∑

j=0j 6=m,n

dj

(

1 + O(1)βj

)

=dmO(1)

dmO(1) + dnO(

1∆x2

)

+r−1∑

j=0j 6=m,n

dj +r−1∑

j=0j 6=m,n

O(1)βj

≤O(1)

O(1) +O(

1∆x2

)

+ (1− dm − dn)= O(∆x2).

Uma característica interessante desta nova fórmula é que o peso atribuído ao estêncil

68

CAPÍTULO 4. O MÉTODO WENO-Z4.1. UMA NOVA FÓRMULA PARA OS PESOS DO WENO

descontínuo é maior do que o peso de Jiang-Shu para o mesmo estêncil. Isto é intuitivo,

já que o peso de Jiang-Shu para estênceis descontínuos é O(∆x4). Podemos verificar este

fato da seguinte maneira: considere quaisquer estênceis Sm e Sn tais que βm ≥ βn. Então

ωzm

ωzn

=dm

(

1 + τ5βm

)

dn

(

1 + τ5βn

) =dmdn

βm

βn

βm + τ5βn + τ5

=dmdn

(

βn

βm

)2(βm

βn

βm + τ5βn + τ5

)

=

=ω(JS)m

ω(JS)n

(

βm

βn

βm + τ5βn + τ5

)

≥ω(JS)m

ω(JS)n

.

Em particular, isto vale para Sm descontínuo. Isto, a princípio, pode parecer uma desvan-

tagem. Afinal, é perigoso atribuir um peso grande à aproximação do estêncil descontínuo,

pois isto pode levar a oscilações e instabilidade numérica. Parece, porém, que há um limite

em que se pode aumentar a contribuição do estêncil descontínuo sem que isto acarrete

problemas. E os pesos de Jiang-Shu estão abaixo deste limite.

Pelo que observamos, o WENO pode lucrar com este pequeno aumento no peso dos

estênceis descontínuos. Quanto mais peso se dá a eles, menos dissipativo é o esquema.

Infelizmente não existe uma estimativa do quanto se pode aumentar do peso sem que isto

gere instabilidade no sistema.

4.1.1 Generalização para todo r

Olhando como o medidor de suavidade global τ5 foi construído, precisamos agora cons-

truir um medidor de suavidade τn que seja uma generalização do caso n = 5, principal-

mente no que diz respeito à optimalidade em regiões contínuas sem pontos críticos. Para

isto, vamos definir os medidores de suavidade τn como o valor absoluto de combinações

69


dos medidores βk,

τn =

∣

∣

∣

∣

∣

r−1∑

k=0

ckβk

∣

∣

∣

∣

∣

. (4.2)

Queremos fazer uma escolha dos ck de modo que

τn = O(∆xr+2), (4.3)

pois, supondo por ora que,

βk = O(∆x2), (4.4)

temos então que

ωzk =

dk

(

1 + τnβk

)

∑r−1j=0 dj

(

1 + τnβj

)

=dk

(

1 + O(∆xr+2)O(∆x2)

)

∑r−1j=0 dj

(

1 + O(∆xr+2)O(∆x2)

)

=dk (1 +O(∆xr))

(1 +O(∆xr))∑r−1

j=0 dj

=dk (1 +O(∆xr))

(1 +O(∆xr))

= dk (1 +O(∆xr)) (1 +O(∆xr)) = dk +O(∆xr).

que é condição suficiente para a optimalidade.

Para manter o caráter global do medidor de suavidade τn, este por sua vez precisa

conter informações sobre todos os pontos de S0, . . . , Sr−1, ou seja,

c0 6= 0, cr−1 6= 0. (4.5)

Neste momento uma pergunta surge: é sempre possível encontrar coeficientes c0, . . . , cr−1

70

CAPÍTULO 4. O MÉTODO WENO-Z 4.2. RESULTADOS AUXILIARES

tais que τn = O(∆xr+2)? Felizmente, a resposta é sim, e veremos este resultado traduzido

em teorema na seção 4.3. Para chegar a este resultado, precisamos antes de tudo encon-

trar alguns resultados auxiliares que nos guiarão à obtenção dos coeficientes que levam à

optimalidade.

4.2 Resultados auxiliares

Antes de entrarmos em detalhes sobre a convergência do método WENO-Z, como a

obtenção dos coeficientes ck em (4.2) de modo que (4.3) seja satisfeita, torna-se neces-

sário obter mais informações sobre os indicadores de suavidade βk. Relembrando por

conveniência a fórmula dos indicadores de suavidade βk,

βk =r−1∑

l=1

∆x2l−1

∫ xi+1

2

xi− 1

2

(

dl

dxlfk(x)

)2

dx,

vemos que para analisar, por exemplo, que βk = O(∆x2), precisamos saber mais sobre

o comportamento dos polinômios fk que aproximam h no estêncil k. É a partir destes

polinômios é que podemos obter resultados significativos a respeito do método WENO-Z.

4.2.1 Representação dos polinômios fk

Por simplicidade, consideremos a partir de agora xi = 0. A partir da equação anterior

vemos que, para obter qualquer informação precisa a respeito dos indicadores βk, precisa-

mos antes de tudo estudar os polinômios fk. Na verdade queremos obter dois resultados

importantes, resumidos no teorema a seguir:

Teorema 6. Seja ρj =⌊

r−j−12

⌋

e f (j) = f (j)(0).Para todo x ∈[

x− 12, x 1

2

]

, fk pode ser

71

4.2. RESULTADOS AUXILIARES CAPÍTULO 4. O MÉTODO WENO-Z

representado como

fk(x) =∑r−1

j=0ak,jx

j, (4.6)

com os coeficientes ak,j expressos como

(a) (Independência)

ak,j =1

j!

ρj∑

δ=0

φ2δf(j+2δ)∆x2δ +O(∆xr−j). (4.7)

ou como

(b) (Simetria)

ak,j =∞∑

l=0

σk,jlf(j+l)∆xl, (4.8)

onde

σk,j,l = (−1)lσr−1−k,jl (4.9)

No caso particular l ≤ ρj,

σk,j,l =

1j!φl mod (l, 2) = 0

0 caso contrário., (4.10)

Antes de demonstrarmos, é necessário entender o que há por trás deste teorema.

Ele nos diz que fk é uma aproximação de h no estêncil k e seus coeficientes ak,j são

aproximações da j-ésima derivada de f no ponto xi. Além disso, expandindo ak,j, vemos

que os termos de ordem até ρj não dependem do estêncil, pois∑ρj

δ=0 φ2δ1j!f (j+2δ)∆x2δ

não depende de k. Quando δ > ρj, os coeficientes de ordem v passam a depender do

estêncil, mas ainda existe a possibilidade de relacioná-los com o seu estêncil simétrico.

A equações (4.8) e (4.9) traduzem o fato de que, ao expandirmos os coefieicentes ak,j

72


do polinômio fk em potências de ∆x, os coeficientes de ordem par serão iguais aos de

ar−1−k,j do polinômio f r−1−k, enquanto os de ordem ímpar terão sinais opostos. Todas

essas propriedades tornam-se fáceis de se visualizar uma vez que observamos a fórmula

explícita de fk para diversos casos. Por exemplo, para n = 5 e k = 0 vemos que

f 0(x) =∑2

j=0a0,jx

j =

=1

24(fj−2 + 2fj−1 + 23fj)

+12(fj−2 − 4fj−1 + 3fj)

∆xx

+1

2

(fj−2 − 2fj−1 + fj)

∆x2x2

Para n = 7 e k = 0,

f 0(x) =1

24(fj−3 − 4fj−2 + 5fj−1 + 22fj)

+124(−7fj−3 + 33fj−2 − 69fj−1 + 43fj)

∆xx

+1

2

(−fj−3 + 4fj−2 − 5fj−1 + 2fj)

∆x2x2

+1

3!

(fj−3 − 3fj−2 + 3fj−1 − 1fj)

∆x3x3. (4.11)

73


Para n = 9 e k = 0,

f 0(x) =1

1920(−71fj−4 + 364fj−3 − 746fj−2 + 684fj−1 + 1689fj)

+148(−9fj−4 + 50fj−3 − 120fj−2 + 174fj−1 − 95fj)

∆xx

+1

2

18(7fj−4 − 36fj−3 + 74fj−2 − 68fj−1 + 23fj)

∆x2x2

+1

3!

12(−3fj−4 + 14fj−3 − 24fj−2 + 18fj−1 − 5fj)

∆x3x3,

+1

4!

(fj−4 − 4fj−3 + 6fj−2 − 4fj−1 + fj)

∆x4x4. (4.12)

Expandindo os coeficientes ak,j em séries de Taylor em torno de 0, podemos observar

claramente que existe uma simetria entre a0,j e a2,j. Notamos também que como a1,j tem

simetria com ele mesmo, os termos de ordem ímpar se cancelam:

a0,0 = f −1

24f ′′∆x2 +

1

24f ′′′∆x3 −

7

288f (4)∆x4 +

1

96f (5)∆x5 (4.13)

a0,1 = f ′ −1

3f ′′′∆x2 +

1

4f (4)∆x3 −

7

60f (5)∆x4 (4.14)

a0,2 =1

2f ′′ −

1

2f ′′′∆x+

7

24f (4)∆x2 −

1

8f (5)∆x3 (4.15)

a1,0 = f −1

24f ′′∆x2 −

1

288f (4)∆x4 (4.16)

a1,1 = f ′ +1

6f ′′′∆x2 +

1

120f (5)∆x4 (4.17)

a1,2 =1

2f ′′ +

1

24f (4)∆x2 (4.18)

74


a2,0 = f −1

24f ′′∆x2 −

1

24f ′′′∆x3 −

7

288f (4)∆x4 −

1

96f (5)∆x5 (4.19)

a2,1 = f ′ −1

3f ′′′∆x2 −

1

4f (4)∆x3 −

7

60f (5)∆x4 (4.20)

a2,2 =1

2f ′′ +

1

2f ′′′∆x+

7

24f (4)∆x2 +

1

8f (5)∆x3 (4.21)

Com esta observação, podemos agora iniciar a demonstração do teorema, que utilizará

o resultado do seguinte lema:

Lema 1. Seja h definido em (2.2). Então

h(x) =∞∑

δ=0

φ2δf(2δ)(x)∆x2δ, (4.22)

O Lema será provado ao final da demonstração do teorema.

Demonstração (a). Lembramos também que fk é uma aproximação de h por um polinô-

mio de grau r, i.e.,

fk(x) = h(x) +O(∆xr). (4.23)

Combinando o resultado do Lema com a equação (4.23) e expandindo f (2δ) em séries de

Taylor em torno de xi, temos

fk(x) =∞∑

δ=0

φ2δ∆x2δf (2δ)(x) +O(∆xr)

=∞∑

δ=0

φ2δ∆x2δ

∞∑

j=0

1

j!f (j+2δ)(0)xj +O(∆xr)

=∞∑

j=0

∞∑

δ=0

φ2δ1

j!f (j+2δ)(0)xj∆x2δ +O(∆xr). (4.24)

75


Como x ∈[

x− 12, x 1

2

]

, x ≤ O(∆x),

φδ

1

j!f (j+2δ)(0)xj∆x2δ ≤ O(∆xj+2δ).

Desta forma, para j+2v > r−1, ou simplesmente para δ > ρj =⌊

r−j−12

⌋

, φδ1j!f (j+2δ)(0)xj∆x2δ ≤

O(∆xr), e portanto

O(∆xr) +∞∑

j=0

∞∑

δ=ρjδ≥0

(

φ2δ1

j!f (j+2δ)(0)xj∆x2δ

)

= O(∆xr). (4.25)

Combinando (4.24) e (4.25),

fk(x) =∞∑

j=0

ρj∑

δ=0

(

φ2δ1


)

+O(∆xr). (4.26)

Observamos também que se j > r − 1, obviamente j + 2δ > r − 1, e com isso,

O(∆xr) +∞∑

j=r

∞∑

δ=0

(

φ2δ1


)

= O(∆xr). (4.27)

De (4.26) e (4.27),

fk(x) =r−1∑

j=0

(

ρj∑

δ=0

φ2δ1


)

+O(∆xr) (4.28)

Colocando xj em evidência no somatório em j,

fk(x) =r−1∑

j=0

(

ρj∑

δ=0

φ2δ1

j!f (j+2δ)(0)∆x2δ +O(∆xr−j)

)

xj,

como queríamos.

76


Demonstração (b). Nosso objetivo é construir aproximações polinomiais de h de grau

(r − 1) nos estênceis Sk e Sr−1−k. Para tal utilizaremos as funções Fk e Gk definidas

abaixo como

Fk(x) =

∫ x

x[k−(r−1)]− 1

2

h(ξ)dξ,

Gk(x) = −

∫ x[(r−1)−k]+ 1

2

x

h(ξ)dξ.

A partir daí obtemos

Fk(xi+ 12) =

∫ xi+1

2

x[k−(r−1)]− 1

2

h(ξ)dξ =i∑

s=k−(r−1)

∫ xs+1

2

xs− 1

2

h(ξ)dξ (4.29)

=i∑

s=k−(r−1)

fs∆x, i ∈ k − (r − 1), ..., k

Gk(xi+ 12) = −

∫ x[(r−1)−k]+ 1

2

xi+1

2

h(ξ)dξ = −

(r−1)−k∑

s=i+1

∫ xs+1

2

xs− 1

2

h(ξ)dξ (4.30)

= −

(r−1)−k∑

s=i+1

fs∆x, i ∈ −k, ..., (r − 1)− k

Seja P Fk (x) o único polinômio de grau r que interpola Fk(x) nos r + 1 pontos xi+ 1

2,

i ∈ k− r, ..., k. Da mesma forma, seja PGk (x) o único polinômio de grau menor ou igual

a r que interpola Gk(x) nos pontos xi+ 12, i ∈ −k− 1, ..., (r− 1)− k. É possível mostrar

que

h(x) = F ′k(x) =

(

P Fk

)′(x) +O(∆xr), x ∈ [x[k−(r−1)]− 1

2, xk+ 1

2],

h(x) = G′k(x) =

(

PGk

)′(x) +O(∆xr), x ∈ [x−k− 1

2, x(r−1)−k+ 1

2],

77


e, desta forma,

fk(x) =(

P Fk

)′(x), x ∈ [x[k−(r−1)]− 1

2, xk+ 1

2],

f (r−1)−k(x) =(

PGk

)′(x), x ∈ [x−k− 1

2, x(r−1)−k+ 1

2],

são as aproximações que queremos. De fato, usando a fórmula de interpolação de Lagrange

adaptada ao nosso contexto,

P Fk (x) =

k∑

i=k−r

Fk(xi+ 12)Lk,i(x),

PGk (x) =

(r−1)−k∑

i=−k−1

Gk(xi+ 12)L(r−1)−k,i(x),

com

Lk,i(x) =k∏

s=k−rs 6=i

x− xs+ 12

xi+ 12− xs+ 1

2

.

Usando (4.29) e (4.30), obtemos

P Fk (x) =

k∑

i=k−r

i∑

s=k−(r−1)

fs∆xLk,i(x)

PGk (x) =

(r−1)−k∑

i=−k−1

(r−1)−k∑

s=i+1

−fs∆xL(r−1)−k,i(x)

=

(r−1)−k∑

i=−k−1

(r−1)−k∑

s=i+1

−fs∆xLk,−i−1(−x)

=k+1∑

i∗=k−(r−1)

(r−1)−k∑

s=−i∗+1

−fs∆xLk,i∗−1(−x)

=k∑

i=k−r

(r−1)−k∑

s=−i

−fs∆xLk,i(−x)

78


Por outro lado,

Lk,i(x) =k∏

s=k−rs 6=i

x− xs+ 12

xi+ 12− xs+ 1

2

=k∏

s=k−rs 6=i

x− (s+ 12)∆x

(i− s)∆x

=

∏k

s=k−r,s 6=i

(

x− (s+ 12)∆x

)

∏k

s=k−r,s 6=i(i− s)∆x=

∑r

j=0 c1k,i,j∆xr−jxj

c2k,i,j∆xr

=r∑

j=0

ck,i,j∆x−jxj

Logo

P Fk (x) =

k∑

i=k−r

i∑

s=k−(r−1)

fs

r∑

j=0

ck,i,j∆x−(j−1)xj

PGk (x) =

k∑

i=k−r

(r−1)−k∑

s=−i

−fs

r∑

j=0

ck,i,j∆x−(j−1)(−x)j

=k∑

i=k−r

i∑

s=k−(r−1)

−f−s

r∑

j=0

ck,i,j∆x−(j−1)(−x)j

Derivando as igualdades acima,

fk(x) =k∑

i=k−r

i∑

s=k−(r−1)

r∑

j=0

fsck,i,j∆x−(j−1)jxj−1

f (r−1)−k(x) =k∑

i=k−r

i∑

s=k−(r−1)

r∑

j=0

−f−sck,i,j∆x−(j−1)(−1)jjxj−1

79


Fazendo Ck,i,j = jck,i,j, usando que (−1)j+1 = (−1)j−1 e reorganizando os índices,

fk(x) =r−1∑

j=0

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

fs

∆xj

xj

f (r−1)−k(x) =r−1∑

j=0

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

f−s

∆xj(−1)j

xj

Desta forma,

ak,j=

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

fs

∆xj

a(r−1)−k,j=

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

f−s

∆xj(−1)j

Expandindo fs em séries de Taylor em torno da origem, e usando que (−1)j = (−1)−j,

ak,j=

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

∑∞δ=0

1δ!f (δ)(0)sδ∆xδ

∆xj,

a(r−1)−k,j=

∑k

s=k−(r−1)

(

∑k

i=s Ck,i,j

)

∑∞δ=0

1δ!f (δ)(0)(−1)δsδ∆xδ

∆xj(−1)j.

Organizando os coeficientes em termos de ∆x,

ak,j=

∞∑

δ=0

k∑

s=k−(r−1)

(

k∑

i=s

Ck,i,j

)

1

δ!f (δ)(0)sδ

∆xδ−j,

a(r−1)−k,j=

∞∑

δ=0

k∑

s=k−(r−1)

(

k∑

i=s

Ck,i,j

)

1

δ!f (δ)(0)(−1)δ−jsδ

∆xδ−j .

80


Fazendo l = δ − j, e constatando através de (4.7) que δ − j ≥ 0,

ak,j=

∞∑

l=0

k∑

s=k−(r−1)

(

k∑

i=s

Ck,i,j

)

sj+l

j + l!f (j+l)(0)

∆xl,

a(r−1)−k,j=

∞∑

l=0

k∑

s=k−(r−1)

(

k∑

i=s

Ck,i,j

)

sj+l

(j + l)!f (j+l)(0)(−1)l

∆xl.

Logo

σk,j,l =k∑

s=k−(r−1)

k∑

i=s

Ck,i,j

sj+l

j + l!

σk,j,l =k∑

s=k−(r−1)

k∑

i=s

Ck,i,j

sj+l

j + l!(−1)l

com σk,jl = (−1)lσr−1−k,jl.

Demonstração do Lema. Usando (2.3), e expandindo h(x±∆x) em x,

f ′(x) = h′(x) +∆x2

4.3!h′′′(x) +

∆x3

16.5!h(5)(x) + ... (4.31)

Integrando ambos os lados,

f(x) = h(x) +∆x2

4.3!h′′(x) +

∆x3

16.5!h(4)(x) + ... (4.32)

Derivando (4.31),

f (2)(x) = h(2)(x) +∆x2

24h(4)(x) + ... (4.33)

f (4)(x) = h(4)(x) +∆x2

24h(6)(x) + ... (4.34)

81


Substituindo h(2) e h(4) em (4.32),

f = h+∆x2

24

(

f (2) −∆x2

24h(4)

)

+∆x5

16.5!h(4) + ... (4.35)

= h+∆x2

24f (2) −

(

1

242−

1

16.5!

)

∆x4

24f (4) + ... (4.36)

= h+1

24f (2)∆x2 −

7

5760f (4)∆x4 + ... (4.37)

e, portanto,

h(x) = f(x)−1

24f (2)∆x2 +

7

5760f (4)∆x4 − ...

Repetindo o processo indefinidamente para ordens maiores de ∆x, obtemos (4.22).

Estes resultados nos dão a partir de agora a possibilidade de analisar a expansão de

fk em potências de ∆x.

4.2.2 Ordem dos indicadores de suavidades βk

Nesta seção veremos um pouco mais a fundo como os indicadores de suavidade βk são

construídos. Para isto, fazemos uso das equações (4.7) e (4.6):

βk =r−1∑

m=1

∆x2m−1

∫ 12∆x

− 12∆x

(

fk(m)(x))2

dx

=r−1∑

m=1

∆x2m−1

∫ 12∆x

− 12∆x

(

r−1−m∑

j1=0

(j1 +m)!

j1!ak,m+j1x

j1

)(

r−1−m∑

j2=0

(j2 +m)!

j2!ak,m+j2x

j2

)

dx

=r−1∑

m=1

∆x2m−1

∫ 12∆x

− 12∆x

r−1−m∑

j1=0

r−1−m∑

j2=0

(j1 +m)!

j1!

(j2 +m)!

j2!ak,m+j1ak,m+j2x

j1+j2dx

=r−1∑

m=1

r−1−m∑

j1=0

r−1−m∑

j2=0

Cm,j1,j2ak,m+j1ak,m+j2∆xj1+j2+2m (4.38)

82


com

Cm,j1,j2 =

(j1+m)!(j2+m)!j1!j2!

2−(j1+j2)

(j1+j2+1)mod (j1 + j2, 2) = 0

0 mod (j1 + j2, 2) = 1, j = 0, . . . , r − 1, l = 0, . . . , r − 1.

Definição 9. Para facilitar a notação do próximo teorema, usaremos a seguinte definição:

E(k,m, j1, j2, l1, l2) = Cm,j1,j2σk,m+j1,l1σk,m+j2,l2 ,

Usando (4.9), obtemos a simetria

E(k,m, j1, j2, l1, l2) = (−1)l1+l2E(r − 1− k,m, j1, j2, l1, l2).

Teorema 7. 4.40βk pode ser escrito como

βk =∞∑

M=2

⌊M2 ⌋∑

j=1

Ak,M,jf(j)f (M−j)∆xM , (4.39)

com

Ak,M,j = (−1)MA(r−1)−k,M,j. (4.40)

Além disso, se j < r e M − j < r, Ak,M,j for independente de k, i.e., então

A0,M,j = A1,M,j... = Ar−1,M,j.

Demonstração. Usando a definição de ak,j em (4.8), obtemos

ak,m+j1ak,m+j2 =∞∑

l1=0

∞∑

l2=0

σk,m+j1,l1σk,m+j2,l2f(m+j1+l1)f (m+j2+l2)∆xl1+l2 .

83


Substituindo na equação (4.38),

βk =r−1∑

m=1

r−1−m∑

j1,j2=0

∞∑

l1,l2=0

E(k,m, j1, j2, l1, l2)f(m+j1+l1)f (m+j2+l2)∆xj1+j2+2m+l1+l2 .

Fixando M1 = m+ j1 + l1 e M2 = m+ j2 + l2, é possível organizar a soma como

βk =∞∑

M1=1

∞∑

M2=1

∑

ΩM1,M2

E(k,m, j1, j2, l1, l2)f(M1)f (M2)∆xM1+M2 ,

onde

ΩM1,M2 = (m, j1, j2, l1, l2) ∈ Ψ1 ∪ . . . ∪Ψr−1 | m+ j1 + l1 = M1 e m+ j2 + l2 = M2,

e

Ψm = 1, ..., r − 1 × 1, ..., r − 1−m2 × N2

Chamando M = M1 +M2,

βk =∞∑

M=2

M−1∑

j=1

∑

Ωj,M−j

E(k,m, j1, j2, l1, l2)f(j)f (M−j)∆xM ,

Como f (j)f (M−j) = f (M−j)f (M−(M−j)), novamente reorganizamos a soma, obtendo

βk =∞∑

M=2

⌊M2 ⌋∑

j=1

∑

Ωj,M−j ∪ ΩM−j,j

E(k,m, j1, j2, l1, l2)f(j)f (M−j)∆xM .

Se

Ak,M,j =∑


E(k,m, j1, j2, l1, l2),

84


nós temos

βk =∞∑

M=2

⌊M2 ⌋∑

j=1

Ak,M,jf(j)f (M−j)∆xM .

Como E(k,m, j1, j2, l1, l2) = (−1)l1+l2E(r − 1 − k,m, j1, j2, l1, l2) e j1 + j2 é sempre par

(caso contrário Cm,j1,j2 = 0), temos que

∑


E(k,m, j1, j2, l1, l2) =∑


(−1)l1+l2E(r − 1− k,m, j1, j2, l1, l2)

=∑


(−1)2m+j1+j2+l1+l2E(r − 1− k,m, j1, j2, l1, l2)

=∑


(−1)ME(r − 1− k,m, j1, j2, l1, l2)

e portanto

Ak,M,j = (−1)MA(r−1)−k,M,j.

Podemos ver também que se j < r e M − j < r, Ak,M,j independe de k. De fato, usando

(4.38) e (4.7),

βk =r−1∑

m=1

r−1−m∑

j1,j2=0

Cm,j1,j2

(ρm+j1∑

δ1=0

ρm+j2∑

δ2=0

φ2δ1φ2δ2∆x2δ1+2δ2

(m+ j1)!(m+ j2)!f (m+j1+2δ1)f (m+j2+2δ2)+

O(∆xr−max(j1,j2)))

∆xj1+j2+2m

A fórmula acima mostra que os coeficientes associados a f (m+j1+2δ1)f (m+j2+2δ2) são inde-

pendentes de k. Desta forma, o maior valor possível para m+j1+2δ1 (e para m+j2+2δ2)

85

4.3. EXISTÊNCIA DE τN CAPÍTULO 4. O MÉTODO WENO-Z

é

m+ j1 + 2ρm+j1 = m+ j1 + 2

⌊

r − (m+ j1)− 1

2

⌋

=

r − 1, se r − (m+ j1) for ímpar,

r − 2, caso contrário.

Sendo assim, se j < r e M − j < r, o coeficiente Ak,M,j independe do estêncil k.

Corolário 3. Se M ≤ r, Ak,M,j é independente de k para todo j ∈ 1, ...,⌊

M2

⌋

.

Teorema 8.

Corolário 4. Seja xi um ponto crítico de ordem ncp = 0. Então βk = O(∆x2).

4.3 Existência de τn

Relembrando a formulação do método WENO-Z o início deste capítulo, vemos que,

para garantir a convergência de ordem n (descrita em (3.11)), é necessário que, dado

uma combinação linear τn =∑r−1

k=0 ckβk, seja satisfeita a τn = O(∆xr+2). Devemos agora

encontrar uma combinação dos ck que satisfaça tal propriedade. Em outras palavras,

expandindo βk em potências de ∆x, queremos que todos os termos de ordem até r+ 1 se

anulem. Fixemos r. Usando (4.39), isto é equivalente a dizer que

r−1∑

k=0

ck

r+1∑

M=2

⌊M2 ⌋∑

j=1

Ak,M,jf(j)f (M−j)∆xM

= 0. (4.41)

86

CAPÍTULO 4. O MÉTODO WENO-Z 4.3. EXISTÊNCIA DE τN

A partir daí,

0 =r−1∑

k=0

ck

r+1∑

M=2

⌊M2 ⌋∑

j=1


=r+1∑

M=2

⌊M2 ⌋∑

j=1

(

r−1∑

k=0

ckAk,M,j

)

f (j)f (M−j)∆xM .

Uma condição suficiente para que a equação anterior seja válida é∑r−1

k=0 ckAk,M,j = 0 para

M = 2, ..., r + 1 e j = 1, ...,⌊

M2

⌋

, então a equação acima é válida. Em outras palavras,

(

A0,M,j A1,M,j · · · Ar−1,M,j

)

.

c0...

cr−1

= 0,M = 2, ..., r + 1,

j = 1, ...,⌊

M2

⌋

.(4.42)

Fixemos M ≥ 2 (note que podemos definir M > r + 1; isto será importante na próxima

seção). Se para todo j a equação anterior for verdadeira, podemos transformar o sistema

de equações em uma multiplicação matricial, com o auxílio da matriz

Br(M) =

A0,(M−1),1 A1,(M−1),1 · · · Ar−1,(M−1),1

A0,(M−1),1 A1,(M−1),1 · · · Ar−1,(M−1),1

......

. . ....

A0,(M−1),⌊M2 ⌋

A1,(M−1),⌊M2 ⌋

· · · Ar−1,(M−1),⌊M

2 ⌋

e, portanto, reescrever (4.42) como

Br(M).

c0...

cr−1

=

0

...

0

, M = 2, ..., r + 1,

87


Juntando todas as matrizes Br(M) em bloco, e construindo

Ar(M) =

Br(3)

...

Br(M)

,

Podemos novamente reescrever (4.42) como

Ar(r + 2).

c0...

cr−1

=

0

...

0

.

Em outras palavras, a equação acima é válida se (c0, ..., cr−1)T for um elemento do

núcleo de Ar(r+2). Na verdade apenas nos interessa um elemento não-trivial do núcleo,

caso contrário τn = 0. Por construção, como os elementos de Ar(r+2) não dependem de

f e suas derivadas, os valores ck também não dependerão.

Observação 6. Uma condição necessária e suficiente para garantir a existência de um

elemento não-trivial da base de Ar(r+2) é que seu posto seja estritamente menor do que

r. De fato, veremos que isto sempre será satisfeito.

Observação 7. Por simplicidade, escreveremos Ar(M) = A(M) e Br(M) = B(M); no

entanto, é importante lembrar que a construção destas matrizes dependem intrinsicamente

da ordem do método WENO-Z utilizado, ou seja, do valor de r. Além disso,

Vamos a alguns exemplos:

Exemplo 2. n = 5(r = 3)

Vimos anteriormente em (3.8) os valores de βk para r = 3. Escolhemos ck da seguinte

88


forma:

c0 = 1, c1 = 0, c2 = −1,

implicando em

τ5 = |β0 − β2| , (4.43)

ou simplesmente

τ5 =

∣

∣

∣

∣

13

3f ′′f ′′′ − f ′f (4)

∣

∣

∣

∣

∆x5 +O(∆x6).

Vemos claramente que a equação acima satisfaz (4.3) e (4.5).

Exemplo 3. n = 7(r = 4)

Em outro exemplo, para r = 4, obtemos os seguintes valores de βk:

β0 = (f ′∆x)2 +13

12f ′′2∆x4 −

1

2f ′f (4)∆x5 +O(∆x6),

β1 = (f ′∆x)2 +13

12f ′′2∆x4 +

1

6f ′f (4)∆x5 +O(∆x6),

β2 = (f ′∆x)2 +13

12f ′′2∆x4 −

1

6f ′f (4)∆x5 +O(∆x6),

β3 = (f ′∆x)2 +13

12f ′′2∆x4 +

1

2f ′f (4)∆x5 +O(∆x6).

Daí obtemos um bom candidato para τ7, que é

τ7 = |β0 − β1 − β2 + β3| , (4.45)

e, consequentemente,

τ7 = O(∆x6).

Note que essa não é a única escolha que satisfaz a equação acima; podemos também

89


escolher

τ7 = |2β0 − 3β2 + β3| .

Poderíamos também escolher

τ7 =∣

∣β70 + β7

1 − 2β72

∣

∣ ,

mas violaríamos a condição (4.5).

Uma análise mais detalhada das possíveis escolhas para τ7 pode ser feita através da

matriz A(7) (com as linhas nulas omitidas) que, neste caso, tem a forma

1 1 1 1

1312

1312

1312

1312

−12

16

−16

12

.

c0

c2

c3

c4

=

0

0

0

0

.

Podemos encontrar uma base ρ1, ρ2 para o seu núcleo, que é dada por

2

0

−3

1

,

1

1

−2

0

.

Qualquer (c0, ..., cr−1)T que seja uma combinação linear destes elementos satisfaz τn =

O(∆xr+2); resta apenas encontrarmos uma combinação que também satisfaça c0 6= 0 e

cr−1 6= 0. Neste caso, uma combinação que atende aos requisitos é (c0, ..., cr−1)T = ρ1−ρ2.

Exemplo 4. n = 9(r = 5)

Repetindo o mesmo processo,

90


β0 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2 −

2

5f ′f (5)

)

∆x6 +O(∆x7),

β1 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2 +

1

10f ′f (5)

)

∆x6 +O(∆x7),

β2 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2 −

1

15f ′f (5)

)

∆x6 +O(∆x7),

β3 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2 +

1

10f ′f (5)

)

∆x6 +O(∆x7),

β4 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2 −

2

5f ′f (5)

)

∆x6 +O(∆x7).

Neste caso é bem fácil de ver diretamente um candidato para τ9,

τ9 = |β0 − β4| , (4.47)

mas como nem sempre é possível, é interessante olharmos para o núcleo de A(9):

0

−1

0

1

0

,

−1

0

0

0

1

,

−13

−23

1

0

0

.

Observamos que a dimensão do núcleo cresce à medida que r aumenta. De fato, é

possível observar que para cada r, a dimensão do núcleo de A(r + 2) é r − 2, e que,

portanto, sempre existirá uma combinação que satisfaça τn = O(∆xr+2). Este fato será

provado no teorema que virá a seguir, mas antes vamos observar mais um caso:

Exemplo 5. n = 11(r = 6)

91


β0 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 −1

3f ′f (6)∆x7 +O(∆x8),

β1 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 +1

15f ′f (6)∆x7 +O(∆x8),

β2 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 −1

30f ′f (6)∆x7 +O(∆x8),

β3 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 +1

30f ′f (6)∆x7 +O(∆x8),

β4 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 −1

15f ′f (6)∆x7 +O(∆x8),

β5 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

360f ′′f (4) +

781

720f ′′′2

)

∆x6 +1

3f ′f (6)∆x7 +O(∆x8),

e um candidato para τ11 é

τ11 = |β0 − β1 − β4 + β5| . (4.49)

Olhando mais atentamente para (4.43) e (4.47), vemos que τ5 e τ9 têm formas muito

similares. Além disso, olhando para (4.45) e (4.49), verificamos que τ7 e τ11 também

apresentam uma forma parecida. Por trás deste fato reside uma propriedade ainda mais

geral que pode ser resumida no seguinte teorema:

Teorema 9. Seja n = 2r − 1 a ordem da reconstrução WENO. Então o indicador de

suavidade de alta ordem

τn =

|β0 − βr−1| mod (r, 2) = 1,

|β0 − β1 − βr−2 + βr−1| mod (r, 2) = 0

é tal que τn = O(∆xr+2).

92


Demonstração. A ideia principal é usar a simetria geométrica e a antissimetria entre fk

e f r−1−k de forma a obter uma relação entre seus coeficientes ak,j and ar−1−k,j.

Usando a definição de βk em (4.39) e a simetria em (4.40), obtemos

β0 =∞∑

M=2

⌊M2 ⌋∑

j=1

A0,M,jf(j)f (M−j)∆xM ,

β1 =∞∑

M=2

⌊M2 ⌋∑

j=1

A1,M,jf(j)f (M−j)∆xM ,

βr−2 =∞∑

M=2

⌊M2 ⌋∑

j=1

A1,M,jf(j)f (M−j)∆xM(−1)M ,

βr−1 =∞∑

M=2

⌊M2 ⌋∑

j=1

A0,M,jf(j)f (M−j)∆xM(−1)M .

Queremos encontrar uma combinação linear dos βk acima de forma que a soma de

todos os termos de ordem O(∆xk) sejam nulos para k < r + 2. Em outras palavras,

procuramos constantes não-triviais c0, c1, cr−2, cr−1 tais que

c0β0 + c1β1 + cr−2βr−2 + cr−1βr−1 = 0,

para todos os termos de ordem até O(∆xr+1). Com isso, obteremos τ = O(∆xr+2).

Usando o Corolário (3), todos os coeficientes de βk associados aos termos ∆xM ,M < r

são iguais.

• Supondo r ímpar, r+1 é par, e os coeficientes de β0 e βr−1 associados ao termo ∆xr+1

são iguais. Além disso, os coeficientes β0 e βr−1 associados aos termos ∆xM ,M ≤

r + 1 são iguais, e portanto

β0 − βr−1 = O(∆xr+2).

93

4.4. ORDEM MÁXIMA DE τN CAPÍTULO 4. O MÉTODO WENO-Z

• Supondo r par, r+1 é ímpar, e o coeficiente associado ao termo ∆xr+1 de βk+βr−1−k

é nulo (devido à simetria encontrada em (4.40)), e para M < r + 1, o coeficiente

de ∆xM é igual a 2∑⌊M

2 ⌋j=1 Ak,M,jf

(j)f (M−j). Isto significa que os coeficientes de

(β0 + βr−1)− (β1 + βr−2) de ordem ∆xM ,M ≤ r + 1 são nulos, proporcionando

β0 − β1 − βr−2 + βr−1 = O(∆xr+2).

4.4 Ordem máxima de τn

Até este ponto, vimos que é possível encontrar τn = O(∆xr+2) de modo a garantir

que a derivada no método WENO-Z seja aproximada com ordem n, desde que βk =

O(∆x2), i.e., que ncp = 0, com ncp definido em (3.9). No entanto, se ncp = 1, temos

pelo primeiro resultado do Teorema 7, que, na vizinhança do ponto crítico, βk = O(∆x4),

e, consequentemente, ω±k − dk = O(∆xr−2), que não é condição suficiente para garantir

ordem n no método. Além disso, se r > 3 , por exemplo, e ncp = 2, βk = O(∆x6) e

portanto, ω±k − dk = O(∆xr−4). A partir deste fato e relembrando a discussão da seção

4.2.2, vemos que à medida que o valor de ncp aumenta, a ordem de convergência do WENO

(não importando se é o WENO clássico ou WENO-Z) deixará de atender a condição de

optimalidade. Em outras palavras, τn = O(∆xr+2) não é mais uma condição que atende

ω±k − dk = O(∆xr) em regiões com pontos críticos. Para contornar este problema, uma

solução é encontrar τn = O(∆xm), m ≥ r + 2 que compense o aumento da ordem de βk.

Desta forma, o que podemos é, para cada r, encontrar M2r−1 = Mn, o maior valor possível

de m tal que existam coeficientes c0, ..., cr−1 que gerem τn = O(∆xMn) 6= 0.

Da mesma forma que fizemos na seção anterior, podemos reescrever τn = O(∆xm)

94

CAPÍTULO 4. O MÉTODO WENO-Z 4.4. ORDEM MÁXIMA DE τN

comor−1∑

k=0

ck

m−1∑

M=2

⌊M2 ⌋∑

j=1


= 0.

e estabelecer uma condição suficiente para que a equação acima seja válida:

A(m).

c0...

cr−1

=

0

...

0

, (c0, ..., cr−1)T 6= 0

Em outras palavras, o que queremos é encontrar o maior valor possível de m tal que a

matriz A(m) ainda tenha uma base não-trivial, i.e., tal que seu posto seja menor do que

r.

Observamos que à medida que o valor de m aumenta, o valor do posto de A(m),

rank(A(m)), também aumenta. Queremos encontrar uma fórmula para esse aumento e

a partir daí uma fórmula fechada para rank(A(m)).

Nossa estratégia se baseia em encontrar rank(A(m + 1)) a partir de rank(A(m)).

Desta forma, para encontrar A(m+ 1), basta adicionar a A(m) as linhas de B(m+ 1):

A(m+ 1) =

A(m)

B(m+ 1)

.

Observação 8. As linhas que serão adicionadas terão impacto direto no aumento do

posto de A(m + 1). No caso em que a linha formada pelos elementos Ak,M,j depende

de k, iremos supor que, ao tomarmos quaisquer duas linhas com essa propriedade, estas

serão linearmente independentes dois a dois. Não provaremos esta suposição, mas ela

parece ser verdadeira, pois não encontramos nenhum contra-exemplo até agora. De fato,

95


a afirmação é válida para r ≤ 20. Se por acaso existir algum contra-exemplo em ordens

maiores que 20, o que poderá acontecer é que para alguns valores de r, o valor encontrado

de M2r−1 será ligeiramente maior, ou seja, um caso melhor do que esperávamos.

Para facilitar, vamos a um exemplo numérico. Com r = 5 , escrevemos os indicadores

de suavidade βk, k = 0, ..., 4, explicitando os termos até ordem 8:

β0 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−2

5f ′f (5) −

1

360f ′′f (4) +

781

720f ′′′2

)

∆x6

+

(

2

3f ′f (6) −

9

5f ′′f (5)

)

∆x7

+

(

−13

21f ′f (7) +

1235

432f ′′f (6) −

5467

1440f ′′′f (5) +

32803

30240f (4)2

)

∆x8 +O(∆x9),

β1 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

1

10f ′f (5) −

1

360f ′′f (4) +

781

720f ′′′2

)

∆x6

+

(

−1

12f ′f (6) +

11

60f ′′f (5)

)

∆x7

+

(

1

21f ′f (7) −

251

2160f ′′f (6) −

781

1440f ′′′f (5) +

32803

30240f (4)2

)

∆x8 +O(∆x9),

β2 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−1

15f ′f (5) −

1

360f ′′f (4) +

781

720f ′′′2

)

∆x6

+

(

−1

126f ′f (7) −

53

2160f ′′f (6) +

781

1440f ′′′f (5) +

32803

30240f (4)2

)

∆x8 +O(∆x9),

β3 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

1

10f ′f (5) −

1

360f ′′f (4) +

781

720f ′′′2

)

∆x6

+

(

1

12f ′f (6) −

11

60f ′′f (5)

)

∆x7

+

(

1

21f ′f (7) −

251

2160f ′′f (6) −

781

1440f ′′′f (5) +

32803

30240f (4)2

)

∆x8 +O(∆x9),

96


β4 = (f ′∆x)2 +13

12f ′′2∆x4 +

(

−2

5f ′f (5) −

1

360f ′′f (4) +

781

720f ′′′2

)

∆x6

+

(

−2

3f ′f (6) +

9

5f ′′f (5)

)

∆x7

+

(

−13

21f ′f (7) +

1235

432f ′′f (6) −

5467

1440f ′′′f (5) +

32803

30240f (4)2

)

∆x8 +O(∆x9),

Com base nas equações acima, escrevemos as matrizes A(m) (com as linhas de valor nulo

omitidas) como

A(6) =

1 1 1 1 1

1312

1312

1312

1312

1312

,

A(7) =

1 1 1 1 1

1312

1312

1312

1312

1312

−25

110

− 115

110

−25

− 1360

− 1360

− 1360

− 1360

− 1360

781720

781720

781720

781720

781720

,

A(8) =

1 1 1 1 1

1312

1312

1312

1312

1312

−25

110

− 115

110

−25

− 1360

− 1360

− 1360

− 1360

− 1360

781720

781720

781720

781720

781720

23

− 112

0 112

−23

−95

1160

0 −1160

95

,

97


A(9) =

1 1 1 1 1

1312

1312

1312

1312

1312

−25

110

− 115

110

−25

− 1360

− 1360

− 1360

− 1360

− 1360

781720

781720

781720

781720

781720

23

− 112

0 112

−23

−95

1160

0 −1160

95

−1321

121

− 1126

121

−1321

1235432

− 2512160

− 532160

− 2512160

1235432

−54671440

− 7811440

7811440

− 7811440

−54671440

3280330240

3280330240

3280330240

3280330240

3280330240

,

Observamos que o posto de cada uma das matrizes é

rank(A(6)) = 1,

rank(A(7)) = 2,

rank(A(8)) = 4,

rank(A(9)) = 5,

e neste caso, portanto, M9 = 8.

Podemos agora discutir como o posto de A(m) varia à medida que o valor de m

aumenta.

Primeiramente, é fácil ver que rank(A(r + 1)) = 1, pois pelo Corolário 3, A0,M,j =

. . . = Ar−1,M,j para M ≤ r, e, portanto, todas as linhas são linearmente dependentes.

Fixemos agora m > r + 1. Para cada m, queremos ver quantas linhas de B(m)

dependem de k. Estas linhas serão as responsáveis pelo aumento do posto de A(m).

98


Pelo Teorema 7, os elementos que não dependem de k são da forma Ak,M,j, com j < r e

m− j < r. Logo o número de linhas novas que não dependem de k é a cardinalidade do

conjunto

j ∈

1, . . . ,⌊m

2

⌋

| j ≥ r ou m− j ≥ r

,

que é rank(B(m)) = m − r (observe que usamos fortemente a suposição contida na

observação [?]).

Em um primeiro momento poderíamos então imaginar que rank(A(m+1)) = rank(A(m))+

rank(B(m + 1)), mas isto não é verdade. Para entender melhor como o posto aumenta,

precisamos analisar as linhas de A(m).

Se r for ímpar, temos que as linhas são da forma

(

A0,m,j A1,m,j · · · A r−12

,m,j · · · Ar−2,m,j Ar−1,m,j

)

.

Com o auxílio da simetria em (2.1),

(

A0,m,j A1,m,j · · · A r−12

,m,j · · · (−1)mA1,m,j (−1)mA0,m,j

)

.

Considerando o espaço formado por estas linhas, vemos que as linhas contém⌈

r2

⌉

variáveis

livres, A0,m,j , ..., A r−12

,m,j se m for par e⌊

r2

⌋

variáveis livres se m for ímpar, pois neste caso

A r−12

,m,j = −A r−12

,m,j = 0, excluindo esta variável livre. Se r for par, não existe termo

central A r−12

,m,j, e portanto há sempre r2

variáveis livres, independentemente do valor de

m.

Ainda no contexto dos espaços formados pelas linhas, vemos através do termo (−1)m

que o k-ésimo elemento de cada linha pode ser igual ou simétrico ao (r − 1 − k)-ésimo

elemento. Essa dualidade nos perimite construir dois espaços distintos, V + e V −, definidos

99


por

V + =

(

a0 a1 · · · ar−1

)

∈ Rr | ak = ar−1−k, k = 0, ..., r − 1

,

V − =

(

a0 a1 · · · ar−1

)

∈ Rr | ak = −ar−1−k, k = 0, ..., r − 1

,

com dimensões⌈

r2

⌉

e⌊

r2

⌋

respectivamente (note que com r par as dimensões são iguais).

Note que V + ∩ V − = 0. Se separarmos A(m) nestes dois espaços, formando duas

novas matrizes A+(m) e A−(m), que contém apenas linhas de V + e V − respectivamente,

temos que

rank(A(m)) = rank(A+(m)) + rank(A−(m)).

Observação 9. Aproveitemos para observar que rank(A+(r + 1)) = 1 e rank(A−(r +

1)) = 0, uma vez que os elementos Ak,m,j independem do estêncil k (vide Corolário 3) e,

portanto, todas as linhas de A(m) são linearmente dependentes e pertencentes a V +.

Com esta separação podemos agora definir rank(A(m+1)) em termos de rank(A(m)).

Com efeito, adicionando as linhas de B(m + 1), vemos que se m for par, apenas o posto

de A+(m + 1) aumentará, pois B(m + 1) ⊂ V +, enquanto que se m for ímpar, apenas

o posto de A−(m + 1) aumentará, pois B(m + 1) ⊂ V −. Por outro lado, esse aumento

está limitado em⌈

r2

⌉

para rank(A+(m+ 1)) e⌊

r2

⌋

para rank(A−(m+ 1)), por conta da

dimensão dos espaços. Desta forma,

rank(A+(m+ 1)) =

min(

rank(A+(m)) + rank(B(m+ 1)),⌈

r2

⌉)

, m par

rank(A+(m)), m ímpar

rank(A−(m+ 1)) =

rank(A−(m)), m par

min(

rank(A−(m)) + rank(B(m+ 1)),⌊

r2

⌋)

, m ímpar

100


Combinando as duas equações acima,

rank(A(m+1)) =

min(

rank(A+(m)) + rank(B(m+ 1)),⌈

r2

⌉)

+ rank(A−(m)), m par

rank(A+(m)) + min(

rank(A−(m)) + rank(B(m+ 1)),⌊

r2

⌋)

, m ímpar

A partir deste resultado obtemos uma relação de recorrência para rank(A(m)), com

condições iniciais rank(A+(r + 1)) = 1 e rank(A−(r + 1)) = 0. A partir do resultado

acima podemos chegar à seguinte conclusão:

Teorema 10. O posto de A(m+ 1), m ≥ r, possui a seguinte forma fechada:

rank(A(m+ 1)) = max(

ρ+(m+ 1),⌈r

2

⌉)

+max(

ρ−(m+ 1),⌊r

2

⌋)

,

com

ρ+(m+ 1) =

1 +⌈

m−r2

⌉2, r ímpar,

1 +(⌊

m−r2

⌋

+ 1) ⌊

m−r2

⌋

, r par,

ρ−(m+ 1) =

(⌊

m−r2

⌋

+ 1) ⌊

m−r2

⌋

, r ímpar,⌈

m−r2

⌉2, r par.

Demonstração. Usaremos o fato de que rank(A+(r + 1)) = 1, rank(A−(r + 1)) = 0 e

rank(B(m)) = m − r. Para cada caso (r par e r ímpar), veremos como os postos de

A+(r + 1) e A

−(r + 1) aumentarão.

Considere inicialmente r ímpar. Usando as relações de recorrência para rank(A+(m+

1)) e rank(A−(m+ 1)), e lembrando que r + 1 é par,

rank(

A+(r + 2)

)

= min(

1 + 1,⌈r

2

⌉)

,

rank(

A−(r + 2)

)

= rank(

A−(r + 1)

)

= 0.

101


De maneira análoga,

rank(

A+(r + 3)

)

= rank(

A+(r + 2)

)

= min(

1 + 1,⌈r

2

⌉)

,

rank(

A−(r + 3)

)

= min(

0 + 2,⌊r

2

⌋)

.

Usando que min(min(a, b) + c, b)) = min(a+ c, b) para a, b, c ≥ 0,

rank(

A+(r + 4)

)

= min(

min(

1 + 1,⌈r

2

⌉)

+ 3,⌈r

2

⌉)

= min(

1 + 1 + 3,⌈r

2

⌉)

,

rank(

A−(r + 4)

)

= rank(

A−(r + 3)

)

= min(

0 + 2,⌊r

2

⌋)

.

rank(

A+(r + 5)

)

= rank(

A+(r + 4)

)

= min(

1 + 1 + 3,⌈r

2

⌉)

,

rank(

A−(r + 5)

)

= min(

min(

0 + 2,⌊r

2

⌋)

+ 4,⌊r

2

⌋)

= min(

0 + 2 + 4,⌊r

2

⌋)

.

Podemos ver que para k ≥ 1, rank (A±(r + k)) assume a forma rank (A+(r + k)) =

min(

ρ+(r + k),⌈

r2

⌉)

e rank (A−(r + k)) = min(

ρ−(r + k),⌊

r2

⌋)

. Basta-nos agora encon-

trar os valores de ρ±(k). Com efeito, organizando os valores encontrados acima na tabela

abaixo, vemos claramente o padrão que se forma:

r + 1 r + 2 r + 3 r + 4 r + 5 r + 6

ρ+(r + k) 1 1 + 1 1 + 1 1 + (1 + 3) 1 + (1 + 3) 1 + (1 + 3 + 5)

ρ−(r + k) 0 0 0 + 2 0 + 2 0 + (2 + 4) 0 + (2 + 4)

Fica evidente que o valor de ρ±(k) aumenta como uma soma de progressão aritmética a

102


cada vez que k é acrescido em 2. Desta forma,

ρ+(r + k) = 1 +

⌊ k2⌋∑

i=1

(2i− 1) = 1 +

⌊

k

2

⌋2

ρ−(r + k) =

⌊ k−12 ⌋∑

i=1

2i =

(

1 +

⌊

k − 1

2

⌋)⌊

k − 1

2

⌋

Escrevendo r + k = m+ 1, substituindo nas equações acima,

ρ+(m+ 1) = 1 +

⌊

m− r + 1

2

⌋2

= 1 +

⌈

m− r

2

⌉2

ρ−(m+ 1) =

(

1 +

⌊

m− r

2

⌋)⌊

m− r

2

⌋

Logo, obtemos

rank(

A+(m+ 1)

)

= min

(

1 +

⌈

m− r

2

⌉2

,⌈r

2

⌉

)

,

rank(

A−(m+ 1)

)

= min

((⌊

m− r

2

⌋

+ 1

)⌊

m− r

2

⌋

,⌊r

2

⌋

)

,

como desejado. Agora consideremos agora o caso r par. Obtemos de forma análoga

(

A+(r + 2)

)

=(

A+(r + 1)

)

= 1,

(

A−(r + 2)

)

= min(

0 + 1,⌊r

2

⌋)

(

A+(r + 3)

)

= min(

1 + 2,⌈r

2

⌉)

,

(

A−(r + 3)

)

=(

A−(r + 2)

)

= min(

0 + 1,⌊r

2

⌋)

103


(

A+(r + 4)

)

=(

A+(r + 3)

)

= min(

1 + 2,⌈r

2

⌉)

,

(

A−(r + 4)

)

= min(

min(

0 + 1,⌊r

2

⌋)

+ 3,⌊r

2

⌋)

= min(

0 + 1 + 3,⌊r

2

⌋)

,

(

A+(r + 5)

)

= min(

min(

1 + 2,⌈r

2

⌉)

+ 4,⌈r

2

⌉)

,

= min(

1 + 2 + 4,⌈r

2

⌉)

(

A−(r + 5)

)

=(

A−(r + 4)

)

= min(

0 + 1 + 3,⌊r

2

⌋)

.

Construindo novamente a tabela,

r + 1 r + 2 r + 3 r + 4 r + 5 r + 6

ρ+(r + k) 1 1 1 + 2 1 + 2 1 + (2 + 4) 1 + (2 + 4)

ρ−(r + k) 0 0 + 1 0 + 1 0 + (1 + 3) 0 + (1 + 3) 0 + (1 + 3 + 5)

obtemos

ρ+(r + k) = 1 +

⌊ k−12 ⌋∑

i=1

2i = 1 +

(

1 +

⌊

k − 1

2

⌋)⌊

k − 1

2

⌋

ρ−(r + k) =

⌊ k2⌋∑

i=1

(2i− 1) =

⌊

k − 1

2

⌋2

Escrevendo r + k = m+ 1, substituindo nas equações acima,

ρ+(m+ 1) = 1 +

(

1 +

⌊

m− r

2

⌋)⌊

m− r

2

⌋

ρ−(m+ 1) =

⌊

m− r + 1

2

⌋2

=

⌈

m− r

2

⌉2

104


(

A+(m+ 1)

)

= min

(

1 +

(

1 +

⌊

m− r

2

⌋)⌊

m− r

2

⌋

,⌈r

2

⌉

)

,

(

A−(m+ 1)

)

= min

(

⌈

m− r

2

⌉2

,⌊r

2

⌋

)

,

como queremos.

A partir do valor obtido para rank(A(m)) podemos encontrar Mn, que é simplesmente

o maior valor possível de m tal que rank(A(m)) < r.

A tabela a seguir contém, para cada n, o valor de Mn:

r n Mn r n Mn

3 5 5 12 23 17

4 7 7 13 25 18

5 9 8 14 27 19

6 11 9 15 29 21

7 13 11 16 31 22

8 15 12 17 33 23

9 17 13 18 35 24

10 19 15 19 37 25

11 21 16 20 39 27

A partir da tabela acima podemos ver que quanto maior for r, a diferença Mn − r,

aumenta, isto é, à medida que r aumenta, existe ainda a possibilidade do esquema garantir

a optimalidade mesmo com a presença de pontos críticos. Por exemplo, quando r = 7,

temos que se ncp = 1, então βk = O(∆x4) e ω±k − dk = O(∆x13−4) = O(∆x9), que ainda

é suficiente para garantir a optimalidade.

105

4.5. UMA ALTERNATIVA PARA OBTER ORDEM NCAPÍTULO 4. O MÉTODO WENO-Z

Abaixo observamos os valores de τn encontrados até ordem 13:

τ5 = |β0 − β2|

τ7 = |β0 + 3β1 − 3β2 − β3|

τ9 = |β0 + 2β1 − 6β2 + 2β3 + β4|

τ11 = |β0 − 10β2 + 10β3 − β5|

τ13 = |β0 + 36β1 + 135β2 +−135β4 − 36β5 − β6|

Pela complexidade das fórmulas para obter estes valores, acreditamos que não existe

uma fórmula fechada que seja simples.

4.5 Uma alternativa para obter ordem n

Observamos nas seções anteriores que quando ncp ≥ 1, nem sempre atingimos a ordem

máxima. Isso ocorre porque à medida que ncp aumenta, a ordem de βk, (denominada

nesta seção por β) também aumenta, mas a ordem de τn continua sendo igual a Mn. Com

isso, a partir de (4.1), vemos que

ω±k − dk = O(∆xMn−β).

ou seja, Mn − β diminui à medida que ncp aumenta. Se Mn − β < r, a condição (3.11)

não é satisfeita e portanto o método não convergirá. Para evitar que isto aconteça, uma

solução aceitável se torna aumentar a ordem de τn, elevando a uma potência p suficiente

para que pMn − β ≥ r. Este valor de p deve ser escolhido o menor possível, pois quanto

maior o valor de p, mais dissipativo o método se tornará. De fato, seja h tal que exista

uma descontinuidade no estêncil Sh. Desta forma, usando βzk como definido em (4.1), e

106

CAPÍTULO 4. O MÉTODO WENO-Z4.5. UMA ALTERNATIVA PARA OBTER ORDEM N

tomando βz∗k como os medidores de suavidade com a razão τn

βkelevado a p, obtemos

βzk

βzh

=βh

βk

βh + τn + ǫ

βk + τn + ǫ<

βh

βk

(βh)q + (τn)

q + ǫ

(βk)q + (τn)

q + ǫ=

βz∗k

βz∗h

,

ou seja, a importância relativa do estêncil Sh diminui à medida que q aumenta, o que em

outras palavras significa tornar o método mais dissipativo.

O valor de p deverá ser escolhido de acordo com o problema em questão. Por exemplo,

suponha r = 6 e ncp = 5. Temos portanto que β = 8 e Mn = 9. Logo

p ≥r + β

Mn

=14

9,

e com p = 149, obtemos

ω±k − dk = O(∆xpMn−β) = O(∆x6),

implicando em convergência de ordem n = 11,

fi+ 12− fi− 1

2

∆x= f ′(xi) +O(∆x11).

Observação 10. Esta potência p pode ser compreendida de maneira análoga à potência

definida para os pesos do WENO clássico, isto é, enquanto no WENO clássico os pesos

são definidos como

ωk =αk

∑r−1j=0 αj

, αk =dk

(βk + ǫ)p,

no método WENO-Z definimos os pesos como

ωzk =

αk∑r−1

j=0 αj

, αk =dkβzk

, βzk =

(

1 +

(

τnβk + ǫ

)p)

107

4.5. UMA ALTERNATIVA PARA OBTER ORDEM NCAPÍTULO 4. O MÉTODO WENO-Z

108

Capítulo 5

Resultados numéricos

No capítulo anterior discutimos sobre o WENO-Z, um esquema alternativo ao WENO

clássico e mapeado, abordados no capítulo 3. Vimos que, dada uma escolha correta do

indicador de suavidade global τn, o método converge com ordem 2r−1 em regiões contínuas

fora da vizinhança de um ponto crítico. Neste capítulo exibiremos o comportamento geral

do esquema WENO-Z em diversos problemas, e compará-los com o WENO clássico e o

WENO mapeado.

Para mostrar o bom comportamento do WENO-Z, em especial o caráter essencialmente

não-oscilatório, vamos utilizar a equação do transporte

ut + ux = 0,

com a condição inicial dada por

u(x, 0) = u0(x).

Utilizamos N = 200 pontos na malha e

∆t = 8∆x53

109

CAPÍTULO 5. RESULTADOS NUMÉRICOS

para que o Runge-Kutta de terceira ordem funcione efetivamente como um de quinta

ordem (em função de ∆x). Neste exemplo utilizaremos uma função descontínua que é

comumente empregada na literatura para este tipo de teste. Ela é a combinação de uma

gaussiana, uma onda quadrada, uma onda triangular e uma semi-elipse. A expressão para

a condição inicial é dada por

u0(x) =

16[G(x, β, z − δ) + 4G(x, β, z) +G(x, β, z + δ)] , x ∈ [−0.8,−0.6]

1 , x ∈ [−0.4,−0.2]

1− [10(x− 0.1)] , x ∈ [0, 0.2]

16[F (x, α, a− δ) + 4F (x, α, a) + F (x, α, a+ δ_)] , x ∈ [0.4, 0.6]

0 , caso contrário,

G(x, β, z) = e−β(x−z)2 ,

F (x , α, a) =√

max(1− α2(x− a)2, 0),

onde as constantes são z = −0.7, δ = 0.005, β = log 236δ2

, a = 0.5 e α = 10.

Calculamos a solução em T = 8 para o método WENO-Z de ordem 5 com p = 1 e 2

como mostra a figura 5.1:

x-1 -0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1Sol. AnalíticaWENO-Z (p=1)WENO-Z (p=2)

Região 1

x-0.4 -0.35 -0.3 -0.25 -0.2

0.92

0.94

0.96

0.98

1

Sol. AnalíticaWENO-Z (p=1)WENO-Z (p=2)

Figura 5.1: Solução da equação do transporte com t = 8 para WENO-Z de ordem 5 comp = 1 e p = 2

110

CAPÍTULO 5. RESULTADOS NUMÉRICOS

Observamos que com p = 1 o método é de fato menos dissipativo, aproximando-se mais

da solução analítica. Poderíamos então imaginar que é sempre melhor escolher p = 1, já

que o método é menos dissipativo. No entanto, em alguns casos, como veremos a seguir,

a escolha de um p pequeno pode acarretar a presença de oscilações. Veremos também em

exemplos posteriores que em ordens maiores esta característica não-oscilatória persistirá,

desde que seja feita uma escolha correta da potência p e do número N de pontos da malha.

Feito este primeiro teste, vamos agora verificar que o método WENO-Z de fato con-

verge com ordem 2r−1 em regiões sem a presença de pontos críticos. Para tal, inicialmente

obtemos uma aproximação da derivada espacial de u em xi, denotada por U(xi). Sendo

u0(xi) a solução analítica de u no ponto xi, podemos medir o erro em xi como sendo

erro(xi) = U(xi)− u0(xi).

Definimos as normas L1, L2, L∞ do erro como

‖erro‖L1=

2

N

N∑

i=0

|erro(xi)| ;

‖erro‖L2=

√

√

√

√

2

N

N∑

i=0

(erro(xi))2;

‖erro‖L∞

= maxi=0,...,N

|erro(xi)| .

Para dois números diferentes de pontos na malha, N1 < N2, é possível obter uma estima-

tiva da ordem do esquema da seguinte maneira: sejam erro1 e erro2 o vetor de erros das

aproximações de um determinado método, usando-se N1 e N2 pontos, respectivamente.

Então a estimativa é dada por

log∥

∥

∥

erro1erro2

∥

∥

∥

log N2

N1

.

111

5.1. ANÁLISE DOS PESOS WENO CAPÍTULO 5. RESULTADOS NUMÉRICOS

No nosso exemplo vamos usar a função exponencial

f(x) = ex

restrita ao domínio [−1, 1] e calcular sua derivada. Neste caso o erro é medido apenas por

erro(xi) = U(xi)− exi .

As tabelas 5.1 a 5.4 mostram a convergência do método WENO-Z para ordens 5, 7, 9 e

11. Podemos ver que os dois esquemas tendem a convergir com ordem 2r − 1.

Observação 11. Utilizamos um número pequeno de pontos para calcular a convergência

do WENO de ordens 9 e 11 devido ao rápido decaimento do erro. Se aumentássemos o

número de pontos da malha, perderíamos a informação da ordem de convergência, uma

vez que o erro se confundiria com o erro da máquina.

5.1 Análise dos pesos WENO

Na seção anterior observamos o caráter não-oscilatório do WENO-Z de ordem 5 através

da função teste, bem como que, em regiões sem pontos críticos, o método converge com

ordem 2r − 1. Precisamos agora comparar estes resultados com os esquemas WENO

clássico (WENO-JS) e WENO mapeado (WENO-M), de modo a verificar em quais casos

o método WENO-Z se torna a melhor opção. Para isso, torna-se necessário fazer uma

análise mais minuciosa dos pesos ωk, primeiramente definidos na seção 3.2.

112

CAPÍTULO 5. RESULTADOS NUMÉRICOS 5.1. ANÁLISE DOS PESOS WENO

N L∞ ordem L2 ordem110 0.21×10−5 0.34×10−6

140 0.55×10−6 5.558 0.80×10−7 5.737170 0.22×10−6 4.828 0.25×10−7 5.655200 0.97×10−7 4.878 0.97×10−8 5.568230 0.49×10−7 4.910 0.45×10−8 5.284

Tabela 5.1: Ordem de convergência do esquema WENO-Z de ordem 5

N L∞ ordem L2 ordem30 0.28×10−10 0.18×10−10

35 0.83×10−11 7.994 0.61×10−11 7.01940 0.31×10−11 7.328 0.24×10−11 7.00745 0.13×10−11 7.521 0.10×10−11 7.00350 0.60×10−12 7.126 0.50×10−12 7.001


N L∞ ordem L2 ordem14 0.38×10−10 0.32×10−10

15 0.20×10−10 8.953 0.17×10−10 8.98516 0.10×10−10 9.017 0.90×10−11 8.98617 0.59×10−11 8.960 0.50×10−11 8.98818 0.34×10−11 9.020 0.29×10−11 8.990


N L∞ ordem L2 ordem14 0.20×10−12 0.17×10−12

15 0.91×10−13 10.92 0.77×10−13 10.9916 0.43×10−13 10.86 0.36×10−13 10.9917 0.22×10−13 10.36 0.18×10−13 10.9518 0.12×10−13 10.83 0.91×10−14 11.04


113


Consideramos a seguinte função:

g(x) =

− sin(πx)− x3

2, se − 1 < x < 0

1− sin(πx)− x3

2, se 0 ≤ x ≤ 1

,

g(x) = g(x− 2), ∀x ∈ R.

S0

S1

S2

Figura 5.2: Função g(x) no intervalo [-1,1]. Os pesos ω0, ω1 e ω2 para a reconstrução dafunção na origem são calculados utilizando os estênceis S0, S1 e S2, respectivamente.

Notamos que g possui uma descontinuidade em x = 0, o que torna esta função inte-

ressante de ser analisada em termos do comportamento dos pesos.

Calculamos os pesos para obter a reconstrução de g (x) utilizando-se o WENO-JS,

WENO-M e WENO-Z, todos de ordem 5, utilizando p = 2. Os resultados podem ser

vistos na figura 5.3.

Podemos reparar que os pesos da parte descontínua (correspondendo aos seis pontos

próximos de x = 0) são ligeiramente maiores no WENO-M do que no WENO-JS e WENO-

Z.

Este resultado pode parecer contraditório para leitor, se comparado com o resultado

obtido em [13], no qual o método WENO-Z aparece com os maiores pesos. Na verdade

114


−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(a) WENO-JS

−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(b) WENO-M

−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(c) WENO-Z

Figura 5.3: Pesos dos esquemas WENO de ordem 5 para o cálculo da função g(x) comp = 2

o que acontece é que a comparação feita foi com o WENO clássico e o WENO mapeado

com p = 2, e WENO-Z com p = 1. Nestas condições o WENO-Z possui pesos maiores

- e consequentemente é menos disispativo. Com efeito, calculando os pesos com p = 1,

e checando o resultado na figura 5.4, vemos que nestas condições os pesos aumentam

significamente.

−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(a) WENO-JS

−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(b) WENO-M

−0.002 −0.001 0 0.001 0.00210

−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

x

ω0

ω1

ω2

(c) WENO-Z


Com este resultado, vemos na verdade que o WENO-M com p = 1 é menos dissipativo

do que o WENO-Z. O método WENO-JS ainda continua com os pesos menores, o que

torna o método mais dissipativo. Confirmamos este fato através das tabelas 5.5 e 5.6.

A tabela nos mostra que em x = −0.001 o estêncil S2 já detectou a descontinuidade

em x = 0, por ser o estêncil mais à direita; por essa razão é que o peso ω2 é bem menor do

que os outros, já que o indicador de suavidade β2 aumenta com a descontinuidade. Além

115


disso em x = 0 observamos que o estêncil S2 já não detém informações a respeito do ponto

de descontinuidade, ao passo que os demais estênceis S0 e S1 detêm. Desta forma seus

pesos diminuem consideravelmente, tornando o peso ω2 maior (na verdade ω2 aumenta

pois é feita uma normalização, aumentando o peso relativo dado a ele).

Podemos também comparar os seis métodos de uma vez (WENO-JS, WENO-M e

WENO-Z para p = 1 e p = 2). A figura 5.5 confirma que a potência p tem grande

influência na dissipação da solução; basta observar que o WENO-JS com p = 1 é menos

dissipativo que WENO-M com p = 2.

x

u

-0.4 -0.35 -0.3 -0.25 -0.2

0.97

0.98

0.99

1

Sol. AnalíticaWENO-M (p=1)WENO-Z (p=1)WENO-JS (p=1)WENO-M (p=2)WENO-Z (p=2)WENO-JS (p=2)

Figura 5.5: Função g(x) no intervalo [-1,1]. Os pesos ω0, ω1 e ω2 para a reconstrução dafunção na origem são calculados utilizando os estênceis S0, S1 e S2, respectivamente.

5.1.1 Ordens mais altas

Podemos estender a discussão da seção anterior para ordens superiores; por exemplo,

para ordem 7, verificamos o seguinte comportamento, apresentado nas figuras 5.6 e 5.7, e

nas tabelas 5.7 e 5.8:

Pelos resultados, observamos novamente que o método WENO-M é o menos dissipativo

116


x = −0.001 x = 0WENO-JS WENO-M WENO-Z WENO-JS WENO-M WENO-Z

ω0 0.1427 0.1271 0.1427 6.5089×10−5 7.1158×10−4 1.3010×10−4

ω1 0.8570 0.8717 0.8569 9.7229×10−4 0.0026 0.0014ω2 2.0045×10−4 0.0011 4.0073×10−4 0.9990 0.9967 0.9985

Tabela 5.5: Medida dos pesos do WENO de ordem 5 com p = 1


ω0 0.1427 0.1272 0.1427 1.2736×10−8 1.4009×10−7 2.5477×10−8

ω1 0.8573 0.8278 0.8573 4.7365×10−7 1.2630×10−6 5.5010×10−7

ω2 9.3794×10−8 5.1386×10−7 1.8767×10−7 0.9999 0.9999 0.9999



ω0 0.0322 0.0322 0.0322 0.9944 0.9789 0.9947ω1 0.3870 0.3867 0.3871 0.0034 0.0129 0.0209ω2 0.5806 0.5807 0.5806 0.0019 0.0055 0.0283ω3 3.6323×10−5 3.9810×10−4 7.2633×10−5 2.9057×10−4 0.0028 0.0062



ω0 0.0322 0.0322 0.0322 0.9999 0.9999 0.9999ω1 0.3870 0.3868 0.3871 9.5771×10−7 3.7510×10−6 3.0159×10−5

ω2 0.5806 0.5810 0.5806 1.9902×10−7 5.8602×10−7 4.4002×10−6

ω3 1.0226×10−8 1.1237×10−7 2.0457×10−8 2.1343×10−8 2.0810×10−7 9.7553×10−6


117


−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

JS (

p=1)

ω0

ω1

ω2

ω3

(a) WENO-JS

−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

M (

p=1)

ω0

ω1

ω2

ω3

(b) WENO-M

−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

Z (

p=1)

ω0

ω1

ω2

ω3

(c) WENO-Z


−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

JS (

p=2)

ω0

ω1

ω2

ω3

(a) WENO-JS

−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

M (

p=2)

ω0

ω1

ω2

ω3

(b) WENO-M

−0.003 −0.002 −0.001 0 0.001 0.002 0.00310

−10

10−9

10−8

10−7

10−6

10−5

10−4

10−3

10−2

10−1

100

xi

Pes

os W

EN

O−

Z (

p=2)

ω0

ω1

ω2

ω3

(c) WENO-Z


dos três, e que a dissipação em cada método pode ser controlada pela potência p utilizada.

Apesar de termos utilizado até então p = 1 e p = 2, estes valores não precisam ser

inteiros. A única questão é que o uso de p não-inteiro torna o método computacionalmente

mais lento, como veremos mais adiante. Veremos também que, apesar do WENO-M ser

o menos dissipativo, o WENO-Z possui um tempo de execução menor que o WENO

mapeado.

Para ilustrar os resultados acima, vamos usar a função teste u0, definida no início deste

capítulo, e comparar os diferentes métodos. Primeiramente, vamos utilizar o WENO de

ordem 7, com N = 200 e p = 2.

Constatamos através do zoom dado na região 1 da figura 5.8 que de fato o WENO-M é

o esquema menos dissipativo dentre os três. Poderíamos diminiuir a potência p de modo

118


x-1 -0.5 0 0.5 1

0

0.2

0.4

0.6

0.8

1Sol. AnalíticaWENO-JSWENO-MWENO-Z

Região 1

x-0.4 -0.35 -0.3 -0.25 -0.2

0.92

0.94

0.96

0.98

1

1.02Sol. AnalíticaWENO-JSWENO-MWENO-Z

Figura 5.8: Solução da equação do transporte com t = 8 para os métodos WENO deordem 7 com p = 2

a tornar o esquema menos dissipativo e, portanto, mais próximo da solução analítica.

Vejamos o que acontece se ao invés de p = 2, utilizarmos p = 1:

x-0.4 -0.35 -0.3 -0.25 -0.2

0.95

1

1.05Sol. AnalíticaWENO-JSWENO-MWENO-Z

Figura 5.9: Zoom na onda retangular da solução da equação do transporte com t = 8para WENO-Z de ordem 7 com p = 1

Ao diminuirmos o valor de p, a solução obtida pelo WENO mapeado apresentou os-

cilações indesejadas. Isto ocorre uma vez que os pesos no estêncil com descontinuidade

não foram pequenos o suficiente, gerando estas oscilações, conhecidas como o fenômeno

119


de Gibbs. Podemos, por outro lado, utilizar uma potência p intermediária.

x-0.36 -0.35 -0.34 -0.33 -0.32 -0.31

0.99

0.995

1

Sol. AnalíticaWENO-M (p=1)WENO-M (p=1.2)WENO-M (p=1.3)WENO-M (p=1.5)WENO-M (p=2)

Figura 5.10: Zoom na onda retangular da solução da equação do transporte em t = 8 doWENO mapeado para diversos valores de p

Neste problema especificamente as oscilações não aparecem a partir de p = 1.3. Este

valor depende do problema e do esquema em questão, da ordem do método e quantos pon-

tos estão sendo utilizados na malha; não faz sentido procurar uma fórmula para encontrar

o melhor p. No entanto é fácil observar que, como o WENO mapeado é menos dissipativo

que o WENO clássico e o WENO-Z, se a potência p escolhida não gerar oscilações no

WENO-M, então não gerará oscilações nos outros dois esquemas.

A primeira impressão que temos é de que então é sempre interessante diminuir o valor

de p ao máximo possível sem gerar oscilações. Isto não é verdade. Utilizar potências

não-inteiras demandam um tempo computacional maior e, dependendo do tamanho do

problema, esta diferença de tempo pode ser crucial para a sua resolução.

120

CAPÍTULO 5. RESULTADOS NUMÉRICOS 5.2. TEMPO DE EXECUÇÃO

p = 1 p = 1.3 p = 1.5 p = 2WENO-JS 33,21s 36,43s 35,40s 34,68sWENO-M 41,66s 42,40s 44,91s 40,21sWENO-Z 31,59s 38,44s 41,94s 34,30s

Tabela 5.9: Tempo de execução da solução do transporte até t = 8 para os esquemasWENO com diferentes valores de p

5.2 Tempo de execução

Discutimos na seção anterior que em determinados casos torna-se interessante utilizar

valores de p não-inteiros, mas que isto demanda um tempo computacional maior. Vamos

confirmar este fato, utilizando para isto o mesmo exemplo da solução do transporte com

a função teste para os esquemas WENO de ordem 7, com N = 2000:

Os três métodos apresentaram um tempo de execução maior com a utilização de p não-

inteiro, mas é notório que o WENO-Z é o que apresentou a maior diferença de resultado

entre as potências. Observamos também uma vantagem do WENO-Z em relação ao

WENO-M em termos computacionais, chegando a ser em média 20% mais eficiente com

p inteiro, e 5% com p não-inteiro. Esta vantagem se deve ao fato do mapeamento feito

no WENO-M ter um alto custo computacional - basta notar que quase não há diferença

entre o WENO clássico e o WENO-Z, especialmente com p inteiro.

5.3 Equações de Euler

Nesta seção, implementaremos o sistema de equações de Euler

ρ

ρv

E

t

+

ρv

ρv2 + p

v(E + p)

x

= 0

121

5.3. EQUAÇÕES DE EULER CAPÍTULO 5. RESULTADOS NUMÉRICOS

com os três esquemas WENO estudados. Neste sistema a variável ρ representa a densidade

de massa, v é a velocidade (ρv é o momento linear), E é a energia total e p é a pressão.

A equação de estado, que relaciona a pressão com as demais variáveis, para os próximos

problemas será dada por

E =p

γ − 1+

1

2ρv2,

com γ sendo uma constante que depende do material do fluido (nestes problemas utiliza-

remos γ = 1, 4). As equações de Euler constituem um sistema não-linear, o que as torna

difíceis de se resolver analiticamente e até numericamente. Além disto, estas equações

são bastante utilizadas na modelagem de fluidos em alta velocidade, onde os termos dis-

sipativos das equações de Navier-Stokes não influenciam muito no resultado. Estas duas

características das equações de Euler fazem delas um teste importante para qualquer

esquema numérico de equações de leis de conservação.

5.3.1 Comparação entre os esquemas WENO

Os próximos testes serão comparações entre os métodos WENO vistos até agora,

usando diferentes potências p. Como em geral não temos solução analítica para estes

problemas, utilizaremos como comparativo o método WENO-M com um número grande

de pontos.

122

CAPÍTULO 5. RESULTADOS NUMÉRICOS 5.3. EQUAÇÕES DE EULER

Interação entre duas blast waves

Este é um problema considerdo difícil na literatura, devido à situação de choques fortes

com ondas de rarefação. A condição inicial é dada por

(ρ, v, p) =

(1, 0, 1000) se 0 ≤ x ≤ 0.1,

(1, 0, 1000) se 0.1 ≤ x ≤ 0.9,

(1, 0, 1000) se 0.9 ≤ x ≤ 1,

com condição de fronteira reflexiva em x = 0 e x = 1.

Calculamos as soluções para o WENO de ordem 9 e p = 5, com tempo final t = 0.038

e N = 300 pontos na malha, exibidas na figura 5.11. Temos duas regiões interessantes que

são analisadas nas figuras adiante. Os detalhes da figura 5.11 confirmam mais uma vez

que o WENO-M é menos dissipativo que o WENO-Z, que por sua vez é menos dissipativo

que o WENO-JS.

Observação 12. Utilizamos o valor p = 5 pois a solução com p = 2 explodiu quando

utilizados o WENO mapeado e o WENO-Z. Na verdade o valor escolhido não é o menor

encontrado de modo a gerar soluções estáveis para os métodos; para o WENO-M o valor

p = 3.3 foi suficiente para gerar uma solução estável, enquanto que para o WENO-Z o

valor de p = 3.9 foi encontrado. Ao contrário do que imaginávamos, o valor encontrado

para o WENO-Z foi maior do que o WENO-M. Deduzimos a partir deste resultado que tal

instabilidade não pode ser justificada somente pela pouca dissipação, mas por outros fatores

que ainda não inspecionamos. Observamos posteriormente que em ordem 7 poderíamos

facilmente utilizar p = 2, e que em ordem 11 nenhum p foi encontrado de modo a gerar

soluções estáveis. Uma alternativa para contornar este problema específico ficará para

trabalhos posteriores.

123


x

Rho

-0.2 0 0.2 0.4 0.6 0.8 1 1.2-1

0

1

2

3

4

5

6

7

WENO-JSWENO-MWENO-ZWENO-M (2000 pontos)

Região 2

Região 1

x

Rho

0.65 0.7 0.75 0.8

4

5

6WENO-JSWENO-MWENO-ZWENO-M (2000 pontos)

x

Rho

0.78 0.8 0.82 0.84 0.86 0.88 0.9

0.8

0.9

1

1.1

1.2

1.3WENO-JSWENO-MWENO-ZWENO-M (2000 pontos)

Figura 5.11: Comparação entre os métodos WENO de ordem 9 no problema da interaçãoentre duas blast waves N = 300 pontos na malha

Problema de Riemann de Lax

Um problema de Riemann é um problema em que a condição inicial é dada por um

choque simples, isto é, há um valor xd tal que

u0(x) =

ul, se x ≤ xd,

ur, se x > xd.

124


No problema de Riemman de Lax a condição inicial é dada por

(ρ, v, p) =

(0.445, 0.698, 3.528) se x ≤ 0,

(0.5, 0, 0.571) se x > 0.

Na figura 5.12 observamos um comparativo entre o método WENO-Z (p = 2) para

ordens 5, 7, 9 e 11 (ou r = 3, 4, 5 e 6, respectivamente). Observamos que à medida que a

ordem do método aumenta, começam a surgir oscilações nas regiões com descontinuidade

(vide região 1 da figura 5.12). Tentamos contornar este problema aumentando a dissi-

pação, isto é, aumentando a potência p, e refinando a malha, sem sucesso. Suspeitamos

que seja devido ao fato desta descontinuidade ser uma descontinuidade de contato, mas é

apenas uma suposição inicial. Esta discussão ficará para trabalhos futuros.

Interação entre um choque e uma onda de entropia

Neste problema a condição inicial é dada por

(ρ, v, p) =

(0.445, 0.698, 3.528) se x ≤ −4,

(0.5, 0, 0.571) se x > −4.

A seguir veremos um panorama geral das soluções obtidas pelo WENO mapeado e WENO-

Z de ordem 11 para a variável ρ, ambos com N = 400 pontos na malha e potência p = 2.

Os dois métodos apresentaram bom comportamento, apesar de não detectarem com

precisão os bicos à esquerda e à direita da oscilação destacada na região 1 da figura

5.13. Observamos também na região 2 que os dois esquemas apresentaram oscilações

indesejadas. Tentamos contornar este problema refinando a malha, como mostra a figura

5.14.

Ao refinarmos a malha, vemos que os dois métodos se aproximaram melhor da solução

125


x

Rho

-4 -2 0 2 4 6

0.4

0.6

0.8

1

1.2

1.4WENO-Z (4000 pontos)WENO-Z (r = 3)WENO-Z (r = 4)WENO-Z (r = 5)WENO-Z (r = 6)

Região 1

Região 2

x

Rho

2 2.5 3

1.24

1.26

1.28

1.3

1.32

1.34

1.36WENO-Z (4000 pontos)WENO-Z (r = 3)WENO-Z (r = 4)WENO-Z (r = 5)WENO-Z (r = 6)

x

Rho

3.2 3.25 3.3 3.35 3.4

0.5

0.51

WENO-Z (4000 pontos)WENO-Z (r = 3)WENO-Z (r = 4)WENO-Z (r = 5)WENO-Z (r = 6)

Figura 5.12: Método WENO-Z com diversas ordens de convergência para o problema deLax N = 400 pontos na malha

real. No entanto, a pouca dissipação do WENO mapeado não foi suficiente para evitar

algumas oscilações, como podemos ver nas regiões 2 e 3 da figura 5.14. Já o WENO-Z se

comportou bem nas três regiões apresentadas, sem apresentar nenhum tipo de oscilação

que destoe da solução desejada.

126


x

Rho

-4 -2 0 2 40

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

WENO-M (2000 pontos)WENO-MWENO-Z

Região 2

Região 1

x

Rho

1 1.5 2 2.5

3

3.5

4

4.5

WENO-M (200WENO-MWENO-Z

x

Rho

-2.9 -2.85 -2.8 -2.75 -2.7 -2.65 -2.6 -2.55 -2.5 -2.45

3.65

3.7

3.75

3.8

3.85

WENO-M (2000 pontos)WENO-MWENO-Z

Figura 5.13: Comparação entre os métodos WENO-M e WENO-Z de ordem 11 no pro-blema da interação entre um choque e uma onda de entropia com N = 400 pontos namalha.

5.3.2 Conclusão

Nos testes exibidos nestas seções vimos que os resultados do WENO-Z foram muito

satisfatórios para diversos problemas clássicos. Em termos de dissipação, o WENO ma-

peado se mostrou o menos dissipativo, seguido do WENO-Z e por fim o WENO clássico.

Vimos também que a potência p utilizada no cálculo dos pesos WENO influi diretamente

127


x

Rho

-4 -2 0 2 40

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

WENO-MWENO-ZWENO-M (2000 pontos)

Região 2

Região 1Região 3

x

Rho

2.3 2.35 2.4

3

3.05

3.1

3.15

3.2

3.25

3.3

3.35


x

Rho

2.38 2.4 2.42 2.44 2.46

0.9

0.95

WENO-JSWENO-MWENO-ZWENO-M (2000 pontos)

x

Rho

-2.9 -2.8 -2.7 -2.6 -2.5

3.65

3.7

3.75

3.8

3.85

3.9


Figura 5.14: Comparação entre os métodos WENO-M e WENO-Z de ordem 11 no pro-blema da interação entre um choque e uma onda de entropia com N = 800 pontos namalha.

128


na dissipação da solução, sendo em alguns casos imprescindível que este valor seja maior

que 1, principalmente em problemas de alta ordem. Em alguns exemplos se o valor de

p fosse pequeno, algumas oscilações indesejadas podem aparecer. Estas oscilações geral-

mente podem ser contornadas amentando esta potência p, ou simplesmente aumentando

o número de pontos da malha. Em outros casos o WENO-Z apresentou oscilações que

não se extinguiram nem com o refinamento da malha nem com o aumento da potência

p. Acreditamos que estas oscilações ocorrem pelo fato da descontinuidade em questão ser

uma descontinuidade de contato, mas é apenas uma suposição. Em termos computacio-

nais, o WENO-Z se mostrou mais rápido do que o WENO-M em média 20% quando p é

inteiro, e que essa melhora cai para 5%, caso contrário.

129


130

Referências Bibliográficas

[1] HARTEN, A. et al. Uniformly high-order accurate essentially non-oscillatory schemes

III. Journal of Computacional Physics, [S.l.], v71, p.231-303, 1987. 5, 7, 8, 39

[2] LIU, X,-D.; OSHER, S.; CHAN T. Weighted essentially non-oscillatory schemes.

Journal of Computacional Physics, [S.l.], v115, 1994, p. 200-212.LIU, X,-D.; OSHER,

S.; CHAN T. Weighted essentially non-oscillatory schemes. Journal of Computacional

Physics, [S.l.], v115, 1994, p. 200-212. 5, 39, 49

[3] HENRICK, A. K.; TASLAM, T. D.; POWERS, J. M. Mapped weighted essentially

non-oscillatory schemes: achieving optimal order near critical points. Journal of Com-

putacional Physics, [S.l.], v207, 2005, p. 542-567 5, 8, 40, 56, 57, 60, 65

[4] JIANG, G.-S.; SHU, C.-W. Efficient implementation of weighted ENO schemes. Jour-

nal of Computacional Physics, [S.l.], v126, 1996, p. 200-228. 5, 8, 13, 19, 39, 40, 49,

56, 65

[5] SHU, C.-W. Essentially non-oscillatory and weighted non-oscillatory schemes for hy-

perbolic conservation laws. ICASE Report, Hampton, no. 97-65, 1997. 12, 32, 35,

50

[6] LEVEQUE, R. J. Numerical methods for conservation laws, Birkhäuser, 2 edition,

Basel, 1992, 232 p. 5, 34

131

REFERÊNCIAS BIBLIOGRÁFICAS REFERÊNCIAS BIBLIOGRÁFICAS

[7] SERNA, S.; MARQUINA, A. Power ENO methods: a fifth-order accurate weighted

power ENO method. Journal of Computational Physics, [S.l.], v194, 2004, p. 632-658.

6

[8] SIDDIQI, K.; KIMIA, B. B.; SHU; C.-W. Geometric shockcapturing ENO schemes

for subpixel interpolation. Graphical Models and Image Processing, 59(5), 1997, p;

278-301. 5

[9] HARTEN, A. Adaptative multiresolution schemes for shock computations. Journal

of Computational Physics, [S.l.], v115, 1994, p.319-338. 37

[10] GOTTLIEB, S.; C.-W. Total variation diminishing Runge-Kutta schemes. Mathe-

matics of computation, Menasha, v. 67, no. 221, 1998, p. 73-85. 19

[11] GELB, A.; TADMOR, E. Detection of edges in spectral data. Applied and computa-

tional harmonic analysis, v7, 1999, p. 101-135. 19

[12] KAMI, S.; KURGANOV, A.; PETROVA, G. A smoothness indicator for adaptative

algorithms for hyperbolic systems. Journal of Computational Physics, [S.l.], v178,

2002, p. 323-341. 19

[13] BORGES, R.; CARMONA, M.; COSTA, B.; DON, W. S. An improved weighted

essentially non-oscillatory scheme for hyperbolic conservation laws. Journal of Com-

putational Physics, [S.l.], v227, 2008, p. 3191-3211. 9, 114

132

Documents

Esquemas Essencialmente Não-Oscilatórios de Alta Precisão ... · Matemática Aplicada, Instituto de Matemática da Universidade Federal do Rio de Janeiro (UFRJ), como parte dos