Embed Size (px)
Citation preview
Centro Federal de Educação Tecnológica de
Minas Gerais
Engenharia Elétrica
MÉTODO ESPECTRAL TAU APLICADO À EQUAÇÃO
DE LAPLACE UNIDIMENSIONAL
Renan Campos Segantini
05/02/2015
Centro Federal de Educação Tecnológica de Minas Gerais Departamento de Engenharia Elétrica Avenida Amazonas, 7675 – Nova Gameleira, Belo Horizonte – MG
(31) 3319-6722
Renan Campos Segantini
MÉTODO ESPECTRAL TAU APLICADO À EQUAÇÃO
DE LAPLACE UNIDIMENSIONA L
Trabalho de conclusão de curso submetida à
banca examinadora designada pelo
Colegiado do Departamento de Engenharia
Elétrica do CEFET-MG, como parte dos
requisitos necessários à obtenção do grau de
bacharel em Engenharia Elétrica.
Área de Concentração: Modelagem
Matemática e computacional
Orientador(a): José Geraldo Peixoto de Faria
Co-orientador(a): Márcio Matias Afonso
CEFET-MG
Belo Horizonte
CEFET-MG
2015
Folha de Aprovação a ser anexada
Aos meus pais, Rosália e Stephenson,
e às minhas irmãs Renata e
Verona.
Agradecimentos
Agradeço ao CEFET, pela valiosa oportunidade de receber tamanho conhecimento.
Ao professor orientador, José Geraldo, que me acompanhou no trabalho, pela
paciência, persistência e grande ajuda no estudo do método espectral tau.
Ao professor co-orientador, Márcio Matias, pelo apoio e grande conhecimento
adquirido através de sua vasta experiência.
Ao meu pai Stephenson, minha mãe Rosália, minhas irmãs Renata e Verona, e a Yuki,
que fizeram de tudo para que eu alcançasse meu sonho.
Agradeço meus colegas do CEFET-MG que, assim como eu, travaram batalhas
intermináveis para conseguir alcançar o sonho de ser engenheiro eletricista.
"...O teu trabalho é a oficina
em que podes forjar a tua própria luz."
Emmanuel
6
Resumo
Métodos de integração numérica cumprem papel de importância quando o
assunto é modelagem matemática e computacional. O avanço das tecnologias e, por
consequência, a melhoria da computação, vem permitindo que os métodos numéricos
sejam cada vez mais viáveis em relação à solução analítica de integrais.
Tendo em vista que uma boa parte dos problemas físicos e de engenharia envolvem
muitas variáveis e que, normalmente, tais problemas são descritos por equações
diferenciais parciais de difícil solução analítica, métodos de integração numérica capazes
de resolver tais equações ganham grande importância.
Buscando a solução numérica de EDPs (equações diferenciais parciais) o objetivo
desse trabalho é o estudo do método espectral tau, método esse independente de malha.
Para a implementação do método espectral tau a revisão bibliográfica traz um
estudo a respeito de outras aplicações de métodos espectrais e também a respeito das
funções de base usadas na interpolação. A implementação do método é feita utilizando
uma ferramenta computacional desenvolvida através do software MATLAB®.
Para validar o método escolheu-se a equação de Laplace unidimensional aplicada
ao problema da distribuição de potencial entre dois pontos, que se caracteriza por um
problema de valor de contorno e por se tratar de um problema de fácil resolução analítica.
7
Abstract
Numerical integration methods have an important role when it comes to
mathematical and computational modeling. The advancement of technologies and,
therefore, the improvement in computing are making numerical methods increasingly
viable when it comes of integrals analytical solution.
Considering that many of the problems regarding physics and engineering involve
many variables and that, usually, such problems are described by partial differential
equations with difficult analytical solution, numerical integration methods capable of
solving such equations become of great relevance.
Looking for numerical solution of PDEs (partial differential equations) the goal of
this paper is the study of spectral tau method, a meshless method.
For the implementation of the spectral tau method, the bibliographic review
presents a study about other applications of spectral methods as well as the basic
functions used in the interpolation. The implementation of this method is done with a
computational tool developed within the MATLAB® software.
To validate the method a one-dimensional Laplace equation applied to the problem
of potential distribution between two points was chosen. This problem have an easy
analytical resolution involving boundary conditions.
.
8
Sumário
Resumo ................................................................................................................................................. 6
Abstract ................................................................................................................................................ 7
Sumário ................................................................................................................................................ 8
Lista de Figuras .............................................................................................................................. 10
Lista de Símbolos ........................................................................................................................... 11
Capítulo 1. Introdução ................................................................................................................ 12
1.1. Relevância do tema em investigação ........................................................................................ 12
1.2. Objetivos do trabalho ..................................................................................................................... 13
1.3. Metodologia ........................................................................................................................................ 13
1.4. Organização do trabalho ............................................................................................................... 14
Capítulo 2. Revisão Bibliográfica ............................................................................................ 15
2.1. Métodos dependentes e independentes de malha e suas características .................. 15
2.2. Introdução aos métodos espectrais .......................................................................................... 16
2.3. Algumas funções de base. (Chebyshev, Hermite e Hermite-Gauss) ............................. 18
2.3.1. Polinômios de Chebyshev...................................................................................................................... 18
2.3.2. Polinômios de Hermite ........................................................................................................................... 19
2.3.3. Funções de Hermite-Gauss ................................................................................................................... 21
2.4. Método Espectral Tau..................................................................................................................... 22
Capítulo 3. Representação Matricial de Operadores e de Produtos Internos .......... 25
3.1. Produto interno e matriz métrica .............................................................................................. 25
3.1.1. Matriz métrica dos polinômios de Chebyshev .............................................................................. 25
3.1.2. Matriz métrica dos polinômios de Hermite ................................................................................... 26
3.1.3. Matriz métrica das funções de Hermite-Gauss ............................................................................. 26
3.2. Operador de derivação .................................................................................................................. 27
3.2.1. Operador de derivação de Chebyshev .............................................................................................. 27
3.2.2. Operador de derivação de Hermite ................................................................................................... 29
3.2.3. Operador de derivação de Hermite-Gauss ..................................................................................... 30
3.3. Operador Laplaciano 𝛁𝟐 ................................................................................................................ 31
9
Capítulo 4. Resultados e Validação do Método ................................................................... 33
4.1. Equação de Laplace e o cálculo da distribuição de potencial entre dois pontos. .... 33
4.2. Método espectral tau aplicado à distribuição de potencial entre dois pontos ......... 35
4.2.1. Resultados obtidos com os polinômios de Chebyshev .............................................................. 38
4.2.2. Resultados obtidos com os polinômios de Hermite ................................................................... 38
4.2.3. Resultados obtidos com as funções de Hermite-Gauss ............................................................. 39
4.3. Validação do método ...................................................................................................................... 40
Capítulo 5. Conclusão de Perspectivas ................................................................................... 41
Referências Bibliográficas ......................................................................................................... 42
Anexos ............................................................................................................................................... 44
10
Lista de Figuras
Figura 2-1 - Em (a) o domínio é preenchido por nós (métodos independentes de malha) e em (b) o domínio
é sub-dividido em malhas (FEM). ....................................................................................................................................... 15
Figura 2-2 - Gráfico dos polinômios de Chebyshev do primeiro tipo até a sexta ordem. ...................................... 19
Figura 2-3 - Polinômios de Hermite do tipo físico até sexta ordem ................................................................................ 21
Figura 2-4 - Função de Hermite-Gauss de grau 0, 10 e 25 e σ=1 ...................................................................................... 22
Figura 4-1 - Geometria do problema. Pontos de potencial determinado. .................................................................... 33
Figura 4-2 - Distribuição de potencial entre dois pontos .................................................................................................... 34
Figura 4-3 - Distribuição de potencial entre dois pontos (polinômios de Chebyshev) .......................................... 38
Figura 4-4 - Distribuição de potencial entre dois pontos (polinômios de Hermite) ............................................... 39
Figura 4-5 - Distribuição de potencial entre dois pontos, 𝜎 = 1 (funções de Hermite-Gauss) ........................... 39
Figura 4-6 - Distribuição de potencial entre dois pontos, 𝜎 = 10 (funções de Hermite-Gauss) ........................ 40
11
Lista de Símbolos
∇² - Operador Laplaciano
∑ - somatório
𝑇𝑘 - Polinômio de Chebyshev de grau k
𝐻𝑘 - Polinômio de Hermite de grau k
∅𝑘 – Função de Hermite-Gauss de ordem k
ℱ - Espaço vetorial onde está definido o produto interno de 𝑇𝑘
𝛿𝑘𝑙 - Deltra de Kronecker
⟨𝑎, 𝑏⟩ - produto interno entre 𝑎 𝑒 𝑏
𝐆𝐜 - Matriz métrica polinômios de Chebyshev
𝐆𝐇 - Matriz métrica polinômios de Hermite
𝐆∅ - Matriz métrica Funções de Hermite-Gauss
𝐆−1 - Inversa da matriz 𝐆
∈ - pertence
�̂�- Operador de derivação
�̃�𝑐 - Matriz que representa o operador �̂� para os polinômios de Chebyshev
�̃�𝐻 - Matriz que representa o operador �̂� para os polinômios de Hermite
�̃�∅ - Matriz que representa o operador �̂� para as funções de Hermite-Gauss
u.c - Unidades de comprimento
EDPs - Equações diferenciais parciais
® - Marca registrada
∀ - Para todo
12
Capítulo 1
Introdução
No Capítulo 1 do texto são abordados a relevância do tema investigado, os
objetivos do trabalho, a metodologia utilizada e a organização do texto.
1.1. Relevância do tema em investigação
Uma boa parte dos fenômenos que cercam nosso dia a dia são compostos de muitas
variáveis. Uma simples previsão do tempo, ondas eletromagnéticas, difusão do calor são
alguns desses problemas que, de certo modo, causam impacto na vida das pessoas. O
desenvolvimento do eletromagnetismo por exemplo, trouxe grandes avanços
tecnológicos para a sociedade. Falar ao celular ou ligar a televisão utilizando controle
remoto, aparelhos de GPS, entre outras tecnologias, envolvem eletromagnetismo e, por
consequência, envolvem métodos capazes de solucionar tais problemas.
Problemas que envolvam um grande número de variáveis são geralmente
descritos por equações diferenciais parciais. Por envolver muitas variáveis, a resolução
analítica de EDPs (equações diferenciais parciais) é, normalmente, muito difícil de se
conseguir. Desta maneira, métodos de integração numérica foram estudados e elaborados
durante os anos e são, até hoje, objeto de estudos, com o objetivo de solucionar problemas
que envolvam EDPs.
Dentre os métodos de integração numérica de equações diferenciais parciais,
existem aqueles que dependem da divisão do domínio de integração em malhas (métodos
com malha ou dependentes de malha) e outros que independem de tal divisão (métodos
sem malha, meshless ou meshfree). Entre os métodos sem malha estão os métodos
espectrais (espectral de Galerkin, espectral tau) que, em geral, podem exibir rápida
convergência justamente por não necessitar de uma subdivisão do domínio de integração
em malhas. Rápida convergência significa normalmente, menor custo computacional para
13
resolver um problema quando se elabora ferramentas computacionais para implementar
tais métodos.
Desta maneira, os métodos espectrais parecem ser interessantes para se estudar,
tendo em vista que podem trazer bons resultados, tanto do ponto de vista matemático
quanto computacional, para a solução de equações diferenciais parciais.
1.2. Objetivos do trabalho
O objetivo geral desse trabalho é estudar métodos de integração numérica capazes
de resolver problemas que envolvam equações diferenciais parciais.
O objetivo específico desse trabalho é estudar o método espectral tau e
implementá-lo à equação de Laplace unidimensional, problema com fácil solução analítica
para que a comparação entre a solução aproximada pelo método e a analítica sirva de
validação.
1.3. Metodologia
Para entender onde os métodos de integração numérica entram, aborda-se no
trabalho a definição de métodos com e sem malha e, para ajudar na compreensão do
método espectral tau, aborda-se a filosofia geral dos métodos espectrais, de onde vieram
e em quais problemas foram implementados com sucesso.
Já visando a implementação do método espectral tau, aborda-se o estudo de
algumas funções que servem de base para o método. Depois do estudo dessas funções, é
feito o estudo do método espectral tau em si e a maneira de representar em forma de
matrizes os produtos internos entre as funções de base e também a representação dos
operadores de derivação de cada função e do operador de Laplace.
Apresenta-se, então, o problema escolhido para validar o método, a distribuição de
potencial entre dois pontos, e a sua solução analítica.
Para implementar o método espectral tau, elabora-se uma ferramenta
computacional baseada no software MATLAB®, escolhido por ser um programa que trata
com facilidade problemas que envolvam matrizes.
14
De posse da solução analítica do problema e da ferramenta computacional,
implementa-se o método espectral tau ao problema utilizando as três funções de base
escolhidas. A comparação entre as duas soluções, analítica e aproximada, é utilizada para
validar o método.
1.4. Organização do trabalho
Esse trabalho está organizado em cinco capítulos e um anexo.
No Capítulo 1, encontra-se a introdução do trabalho onde se discute a relevância
do tema investigado, os objetivos do trabalho, a organização e metodologia utilizada.
No Capítulo 2, apresenta-se uma introdução sobre métodos de integração
numérica com e sem malha, apresenta-se os métodos espectrais em geral, origem e
aplicações em outras áreas e suas considerações matemáticas. Ainda, no Capítulo 2, tem-
se a descrição geral do método espectral tau, onde surgiu e em quais problemas foi
implementado. Apresenta-se uma descrição geral sobre as funções de base, os polinômios
de Chebyshev, os polinômios de Hermite e as funções de Hermite-Gauss, usadas para fazer
a interpolação.
No Capítulo 3, aborda-se a representação matricial dos produtos internos das
funções de base, a representação matricial dos operadores de derivação de cada função e
a representação do operador de Laplace para o caso unidimensional.
O Capítulo 4 estão os resultados. Nele constam a solução analítica da equação de
Laplace aplicada ao cálculo da distribuição de potencial entre dois pontos e a
implementação do método espectral tau ao mesmo problema. No Capítulo 4 encontra-se
também a validação do método na comparação entre as soluções analíticas e
aproximadas.
No Capítulo 5 estão as conclusões sobre a implementação do método e as
perspectivas para trabalhos futuros.
No Anexo estão os códigos da ferramenta computacional desenvolvida em
MATLAB® para a implementação do método.
15
Capítulo 2
Revisão Bibliográfica
2.1. Métodos dependentes e independentes de malha e suas características
Dentre os métodos numéricos mais utilizados atualmente para resolver problemas
que envolvam EDP’s (equações diferenciais parciais) é possível citar o FEM (método dos
elementos finitos), um método dependente de malha. Nos métodos dependentes de malha
há a necessidade de se dividir o domínio de integração em subdomínios para aproximar a
solução. Em casos onde se usa o recurso computacional para alcançar a solução, essa
subdivisão do domínio pode exigir alto custo de memória e processamento por haver a
necessidade de se otimizar a distribuição das malhas em casos onde existem quinas e
grandes variações da grandeza de interesse.
Os métodos independentes de malha (ou sem malha), como o próprio nome indica,
não subdivide o domínio de integração em malhas menores mas espalha pontos (nós)
sobre o contorno (fronteira) e no interior do domínio (nuvem de pontos). A Figura 2-1 a
seguir mostra como o domínio é interpretado por cada método:
Figura 2-1 - Em (a) o domínio é preenchido por nós (métodos independentes de malha) e em (b) o domínio é subdividido em malhas (FEM).
Nos métodos dependentes de malha, a aproximação da solução é calculada
baseando-se em interpolações locais em cada elemento. Dessa maneira, a aproximação é
válida somente para tal subdomínio e a continuidade da aproximação se garante quando
há uma conexão pré-determinada entre essas divisões. Já nos métodos independentes de
malha, a aproximação é feita em cada posição do domínio como um todo, de maneira
16
dinâmica, mas sem perder seu caráter local onde cada posição causa influência nas
proximidades. Assim é possível estabelecer novas relações de conectividade entre os
pontos através da distribuição de nós neste (DUARTE; ODEN 1995).
Dentre os métodos independentes de malha pode-se citar, como exemplos:
Método dos elementos difusos (NAYROLES; TOUZOT; VILLON 1992);
Método de Galerkin Livre de Elementos (BELYTSCHKO; LU; GU 1994);
Métodos espectrais (ORSZAG; GOTTLIEB);
2.2. Introdução aos métodos espectrais
Os métodos espectrais, de maneira geral, tiveram suas primeiras aplicações na
solução de EDPs (equações diferenciais parciais), no fim da década de 60, com os
trabalhos de Orszag e Eliassem, Machenhauer e Rasmussen.
Nos trabalhos de Orzag e Gottlieb, consta que os métodos espectrais são aplicáveis
com grande sucesso em mecânica dos fluidos onde a hidrodinâmica é usada para o estudo
de turbulências, previsão do tempo e a dinâmica dos oceanos[1].
O que dá nome ao método (espectral) é a propriedade de sua expansão em série,
em que o erro de truncamento entre a série com um número finito de termos, N, e a função
a ser aproximada decai a zero mais rápido que qualquer potência de 1/N por ser
infinitamente diferenciável ou analítica. Esse tipo de convergência é chamada de
convergência espectral[2].
Nos métodos espectrais, busca-se a aproximação de um conjunto de EDPs
(equações diferenciais parciais) através de uma interpolação usando funções de base
ortogonais. As expansões em série de funções ortogonais caracterizadas pelos métodos
espectrais são infinitamente diferenciáveis e são rapidamente convergentes caso não
existam variações bruscas nas condições do problema (condições de contorno por
exemplo). Esses métodos buscam a minimização dos resíduos resultantes da resolução
das equações diferenciais e que, dependendo do tipo de minimização, podem ser
classificados da seguinte maneira:
a) Métodos pseudo-espectrais
17
b) Método espectral de Galerkin
c) Método espectral tau
Nos métodos pseudo-espectrais (a), obtém-se resíduos nulos em pontos de
colocação previamente escolhidos. No espectral de Galerkin (b), as funções teste são
iguais às funções tentativa e ortogonais ao resíduo da equação diferencial e no método
espectral tau (c), o resíduo é também ortogonal às funções teste, mas essa condição é
imposta para um subconjunto do total das funções teste e os graus de liberdade restantes
são empregados para impor as condições de contorno. Os métodos espectrais são muito
usados como alternativa aos métodos das diferenças finitas e dos elementos finitos para
resolver equações diferenciais. Dentre as vantagens dos métodos espectrais, está a
produção de uma solução global, com rápida convergência[3].
Os métodos espectrais, basicamente, aproximam a solução procurada através de
uma combinação linear, ou interpolação, de outras funções pré-definidas. Para melhor
compreender tal procedimento, imagine que se deseje conhecer a função 𝑢(𝑡), função essa
que é a solução de um sistema de equações diferenciais. Os métodos espectrais vão
calcular uma outra função, 𝑃(𝑡), tal que 𝑃(𝑡) ≃ 𝑢(𝑡). Mas, conforme dito anteriormente, a
função que aproxima a solução do sistema, 𝑃(𝑡), deve ser composta por uma combinação
linear de outras funções, então,
𝑢(𝑡) = 𝑃(𝑡) =∑𝑎𝑖
∞
𝑖=1
𝑓𝑖(𝑡). (2.2)
Os coeficientes 𝑎𝑖 da expansão, são as incógnitas a serem determinadas através do
método espectral escolhido, enquanto as funções 𝑓𝑖(𝑡), são previamente escolhidas como
interpoladores dessa expansão e são as funções de base da aproximação. Para que a
expansão seja precisamente igual a solução procurada, a combinação linear das funções
de base deve ser infinita, mas tal situação não é factível. Então, normalmente, é realizado
o truncamento dessa combinação em um número inteiro 𝑁, e esse truncamento introduz
um resíduo que faz com que 𝑃(𝑡) não seja precisamente igual a 𝑢(𝑡). Ao realizar o
truncamento da série obtém-se:
𝑢(𝑡) ≃ 𝑃(𝑡) =∑𝑎𝑖
𝑁
𝑖=1
𝑓𝑖(𝑡), 𝑖 = 1,2… ,𝑁,
e tal truncamento introduz um resíduo, 𝑟(𝑡), tal que:
18
𝑟(𝑡) = 𝑃(𝑡) − 𝑢(𝑡).
Quando se usa aproximações numéricas para problemas, sempre há busca pela
minimização do resíduo e, para isso, usa-se um alto valor para N. Porém, quanto maior for
o valor de N, maior será o custo computacional para o cálculo da aproximação (quando a
solução o método é aplicado através de ferramenta computacional). Dessa maneira, deve-
se buscar atender os dois aspectos, tanto a precisão do método quanto o baixo custo de
processamento.
2.3. Algumas funções de base. (Chebyshev, Hermite e Hermite -Gauss)
Em 2.3 serão mostradas algumas funções de base, ortogonais, passíveis de
utilização como interpoladores em métodos espectrais. Elas são os polinômios de
Chebyshev, os polinômios de Hermite e as funções de Hermite-Gauss.
2.3.1. Polinômios de Chebyshev
A grande vantagem do uso dos polinômios de Chebyshev como base para a
aproximação é que essas funções são fáceis de ser calculadas dentro do domínio para os
quais são propostos. Os polinômios de Chebyshev são empregados na expansão de
funções não-periódicas.
Os polinômios de Chebyshev de ordem 𝑘, ⧼𝑇𝑘(𝑥)⧽𝑘=0,1,… , são definidos como
𝑇𝑘(𝑥) = cos(𝑘 arccos x) e o produto interno entre eles, definido na base dos polinômios
de Chebyshev é [4]
⟨𝑇𝑘, 𝑇𝑙⟩ = ∫𝑇𝑘(𝑥)𝑇𝑙(𝑥)𝑑𝑥
√1 − 𝑥2= 𝑐𝑘𝛿𝑘𝑙
π
2, (2.3.1)
1
−1
sendo 𝛿𝑘𝑙 o delta de Kronecker definido por
𝛿𝑘𝑙 = {1, 𝑘 = 𝑙0, 𝑘 ≠ 𝑙
e
19
𝑐𝑘 = {2, 𝑘 = 01, 𝑘 ≠ 0
.
Os 7 (sete) primeiros polinômios de Chebyshev do primeiro tipo são:
𝑇0(𝑥) = 1
𝑇1(𝑥) = 𝑥
𝑇2(𝑥) = 2𝑥2 − 1
𝑇3(𝑥) = 4𝑥3 − 3𝑥
𝑇4(𝑥) = 8𝑥4 − 8𝑥2 + 1
𝑇5(𝑥) = 16𝑥5 − 20𝑥3 + 5𝑥
𝑇6(𝑥) = 32𝑥6 − 48𝑥4 + 18𝑥2 − 1
e seus gráficos podem ser vistos na Figura 2-2 [5]:
Figura 2-2 - Gráfico dos polinômios de Chebyshev do primeiro tipo até a sexta ordem. [5]
2.3.2. Polinômios de Hermite
Os polinômios de Hermite são funções de base ortogonal (assim como os de
Chebyshev). Esses polinômios têm grande aplicação em mecânica quântica e são
especialmente utilizados em estudos com osciladores harmônicos unidimensionais. Essas
funções ganharam tal nome em homenagem a Charles Hermite[6].
Existem dois tipos de polinômios de Hermite, os “polinômios de Hermite
probabilísticos” e os “polinômios de Hermite físicos”. O primeiro é definido pela seguinte
relação:
𝐻𝑛(𝑥) = (−1)𝑛𝑒𝑥2
2𝑑𝑛
𝑑𝑥𝑛𝑒−
𝑥2
2 .
20
Já os polinômios de Hermite físicos (utilizados especificamente nesse trabalho) são
definidos pela seguinte relação:
𝐻𝑛(𝑥) = (−1)𝑛𝑒𝑥
2 𝑑𝑛
𝑑𝑥𝑛𝑒−𝑥
2 (2.3.2.1),
onde 𝑛 é o grau do polinômio. [6]
Seja o polinômio de Hermite de grau 𝑛, 𝐻𝑛(𝑥) e outro de grau 𝑚,𝐻𝑚(𝑥), o produto
interno entre eles é definido por:
⟨𝐻𝑚, 𝐻𝑛⟩ = ∫ 𝐻𝑚(𝑥)𝐻𝑛(𝑥)𝑒−𝑥2𝑑𝑥
+∞
−∞
(2.3.2.2),
onde 𝑒−𝑥2 é a função peso do produto interno. A primeira diferença notória entre os
polinômios de Hermite e os polinômios de Chebyshev, é que os de Hermite são definidos
para toda a reta real enquanto os de Chebyshev são definidos para o intervalo de -1 a 1.
Assim como os polinômios de Chebyshev, os polinômios de Hermite formam uma base
ortogonal. Assim, o seu produto interno é
⟨𝐻𝑚, 𝐻𝑛⟩ = {0, 𝑠𝑒 𝑚 ≠ 𝑛
2𝑛𝑛! √𝜋 , 𝑠𝑒 𝑚 = 𝑛 (2.3.2.3).
Os 7 primeiros polinômios de Hermite do tipo físico são
𝐻0(𝑥) = 1
𝐻1(𝑥) = 2𝑥
𝐻2(𝑥) = 4𝑥2 − 2
𝐻3(𝑥) = 8𝑥3 − 12𝑥
𝐻4(𝑥) = 16𝑥4 − 48𝑥2 + 12
𝐻5(𝑥) = 32𝑥5 − 160𝑥3 + 120𝑥
𝐻6(𝑥) = 64𝑥6 − 480𝑥4 + 720𝑥2 − 120
Na Figura 2-3 são apresentados os polinômios acima descritos para um intervalo
entre -1 e 1.
21
Figura 2-3 - Polinômios de Hermite do tipo físico até sexta ordem
2.3.3. Funções de Hermite-Gauss
As funções de Hermite-Gauss são funções construídas a partir dos polinômios de
Hermite. O que acontece em tais funções é a introdução de um contorno gaussiano ao
redor do polinômio de Hermite desejado, cuja largura definida por um fator 𝜎.
As funções de Hermite-Gauss são dadas por [7]
∅𝑛(𝑡) =1
√2𝑛𝑛! 𝜎√𝜋𝑒−
(𝑡/𝜎)2
2 𝐻𝑛 (𝑡
𝜎) (2.3.3.1),
onde 𝐻𝑛(𝑡) são polinômios de Hermite de ordem 𝑛 e 𝜎 é o fator de escala, fator esse que
define a largura do contorno gaussiano.
O produto interno entre funções de Hermite-Gauss é definido por,
⟨∅𝑛(𝑡), ∅𝑚(𝑡)⟩ =1
√2𝑛𝑛!𝜎√𝜋
1
√2𝑚𝑚!𝜎√𝜋∫ 𝐻𝑛 (
𝑡
𝜎)𝐻𝑚 (
𝑡
𝜎) 𝑒−
(𝑡/𝜎)2
2 𝑑𝑡 (2.3.3.2)+∞
−∞,
e esse produto interno é
⟨∅𝑚(𝑡), ∅𝑛(𝑡)⟩ = {0, 𝑠𝑒 𝑛 ≠ 𝑚1, 𝑠𝑒 𝑛 = 𝑚
(2.3.3.3).
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-150
-100
-50
0
50
100
150
200Polinômios de Hermite Fisicos
X
Hn(X
)
n=0
n=1
n=2
n=3
n=4
n=5
n=6
22
Utilizando a Definição (2.3.3.3), sabe-se então que as funções de Hermite-Gauss
definidas em (2.3.3.1) formam uma base ortonormal.
A Figura 2 mostra os valores assumidos por ∅𝑛(𝑡) para 𝑛 = 0, 𝑛 = 10 e 𝑛 = 25
para 𝜎 = 1,
Figura 2-4 - Função de Hermite-Gauss de grau 0, 10 e 25 e σ=1
2.4. Método Espectral Tau
O método tau foi primeiramente proposto por Lanczos [8] e seu uso com os
polinômios de Chebyshev (funções de Chebyshev do primeiro tipo) foi desenvolvido por
Fox [9] sendo aplicado e defendido por Orszag para uma ampla variedade de problemas
[1]. O método tau utiliza uma expansão em série truncada em um conjunto de funções
completas como uma aproximação para a solução de uma equação diferencial ordinária.
O método também pode ser aplicado em equações diferenciais parciais, equações
integrais e equações integro-diferenciais. [1]
Para aplicar o método espectral tau a um problema de valor de contorno, deve-se
incluir à expansão modos correspondentes às condições de contorno. Dessa maneira, os
modos que devem ser incluídos são numericamente iguais à quantidade de condições de
contorno, então, para 𝑁𝑏 condições de contorno, deve-se adicionar 𝑁𝑏 modos[10].
Suponha que, para o intervalo −𝑎 ≤ 𝑥 ≤ 𝑎 deseja-se resolver a seguinte equação
diferencial
𝑑𝑓(𝑥)
𝑑𝑥= 0,
-10 -8 -6 -4 -2 0 2 4 6 8 10-0.6
-0.4
-0.2
0
0.2
0.4
0.6
0.8
Tempo (s)
Herm
Gauss(n
,t)
Funções de Hermite-Gauss para Sigma=1
n=0
n=10
n=25
23
e com a seguinte condição de contorno,
𝑓(−𝑎) = 𝑉1.
Para resolver o sistema de equações diferenciais e algébricas acima utilizando o
método espectral tau, deve-se aproximar a função 𝑓(𝑥) por um somatório com 𝑁 +𝑁𝑏
modos, sendo 𝑁𝑏 os modos adicionados pelas condições de contorno. Dessa maneira,
sabe-se que 𝑁𝑏 = 1 pois o problema envolve uma condição de contorno. Sendo 𝑃(𝑥) a
função que aproxima 𝑓(𝑥), temos
𝑓(𝑥) ≈ 𝑃(𝑥) = ∑ 𝑢𝑘𝑤𝑘(𝑥),
𝑁+𝑁𝑏
𝑘=0
(2.4.1)
onde 𝑤𝑘(𝑥) são as funções de base da interpolação. Expandindo também as condições de
contorno obtém-se
𝑓(−𝑎) ≈ 𝑃(−𝑎) = ∑ 𝑢𝑘𝑤𝑘(−𝑎) = V1
𝑁+𝑁𝑏
𝑘=0
(2.4.2).
De acordo com a equação diferencial, deve-se, então, aplicar o operador de
derivação também à Equação (2.4.1) para validar o sistema imposto. Assim,
𝑑𝑓(𝑥)
𝑑𝑥≈𝑑𝑃(𝑥)
𝑑𝑥=𝑑
𝑑𝑥( ∑ 𝑢𝑘𝑤𝑘(𝑥)
𝑁+𝑁𝑏
𝑘=0
) = 0.
Como os 𝑢𝑘 são constantes (a serem determinadas), é possível passá-las para fora
do operador de derivação e, juntando com a Equação (2.4.2) obtém-se o seguinte sistema:
{
𝑑𝑃(𝑥)
𝑑𝑥= ∑ 𝑢𝑘
𝑁+𝑁𝑏
𝑘=0
𝑑𝑤𝑘(𝑥)
𝑑𝑥= 0
𝑃(−𝑎) = ∑ 𝑢𝑘𝑤𝑘(−𝑎) = V1
𝑁+𝑁𝑏
𝑘=0
Agora, o passo para obter o sistema a ser resolvido, é tomar o produto interno
entre a primeira equação do sistema anterior pela mesma função de base, 𝑤𝑙, com 0 ≤ 𝑙 ≤
𝑁, e substituir a última equação obtida de tal produto, pela equação que trata a condição
de contorno. Obtém-se, então,
24
{
⟨𝑤𝑙 ,
𝑑𝑃(𝑥)
𝑑𝑥⟩ = ⟨𝑤𝑙 , ∑ 𝑢𝑘
𝑁+𝑁𝑏
𝑘=0
𝑑𝑤𝑘(𝑥)
𝑑𝑥⟩ = 0; 0 ≤ 𝑙 ≤ 𝑁
𝑃(−𝑎) = ∑ 𝑢𝑘𝑤𝑘(−𝑎)
𝑁+𝑁𝑏
𝑘=0
= 𝑉1;
e tal sistema é formado, então, por N+1 incógnitas e N+1 equações que pode ser facilmente
resolvido por se tratar de um sistema linear[1].
25
Capítulo 3
Representação Matricial de Operadores e de Produtos Internos
No Capítulo 3 estão as representações matriciais dos produtos internos das
funções de base apresentadas no Capítulo 2, como também os operadores de derivação e
o operador Laplaciano (derivada segunda).
É importante ressalvar que tudo que diz respeito à matriz aparece em negrito.
Assim o termo genérico [𝑨]𝒊,𝑗 é o elemento da matriz A, que está na linha “i” e na coluna
“j”.
3.1. Produto interno e matriz métrica
O produto interno entre quaisquer funções de base pode ser representado
matricialmente. A matriz que representa o produto interno entre tais funções é chamada
matriz métrica. Os produtos internos entre as funções de base descritas no Capítulo 2 são
calculadas de maneiras distintas para cada função, dessa forma, cada uma delas possui
uma matriz métrica diferente.
3.1.1. Matriz métrica dos polinômios de Chebyshev
De acordo com a Equação (2.3.1), o produto interno entre um polinômio de Chebyshev de
grau k, 𝑇𝑘(𝑥), e um de grau l, 𝑇𝑙(𝑥), é dado por,
⟨𝑇𝑘 , 𝑇𝑙⟩ = 𝑐𝑘𝛿𝑘𝑙π
2,
sendo
𝛿𝑘𝑙 = {1, 𝑘 = 𝑙0, 𝑘 ≠ 𝑙
e
26
𝑐𝑘 = {2, 𝑘 = 01, 𝑘 ≠ 0
.
Com isso é possível definir a matriz métrica dos polinômios de Chebyshev ([𝑮𝒄]𝒌,𝒍)
como,
⟨𝑇𝑘 , 𝑇𝑙⟩ = [𝑮𝒄]𝑘,𝑙 (3.1.1)
onde
𝑮𝒄 =
[ 𝜋 0 … 0
0𝜋
2… 0
⋮ ⋮ ⋱ ⋮
0 0 ⋯𝜋
2]
.
3.1.2. Matriz métrica dos polinômios de Hermite
De acordo com a Equação (2.3.2.3), o produto interno entre um polinômio de
Hermite de grau m, 𝐻𝑚(𝑥) e outro de grau n, 𝐻𝑛(𝑥), é dado por,
⟨𝐻𝑚, 𝐻𝑛⟩ = {0, 𝑠𝑒 𝑚 ≠ 𝑛
2𝑛𝑛! √𝜋 , 𝑠𝑒 𝑚 = 𝑛.
Com isso é possível definir a matriz métrica dos polinômios de Hermite ([𝑮𝑯]𝑚𝑛)
como,
⟨𝑇𝑚 , 𝑇𝑛⟩ = [𝑮𝑯]𝒎,𝒏 (3.1.2)
onde,
𝑮𝑯 =
[ 211! √𝜋 0 ⋯ 0
0 222! √𝜋 ⋯ 0⋮ ⋮ ⋱ ⋮0 0 ⋯ 2𝑛𝑛! √𝜋]
.
3.1.3. Matriz métrica das funções de Hermite-Gauss
De acordo com a Equação (2.3.3.3), o produto interno entre um função de Hermite-
Gauss de grau m, ∅𝑚(𝑥) e outra de grau n, ∅𝑛(𝑥), é dado por,
⟨∅𝑚(𝑥), ∅𝑛(𝑥)⟩ = {0, 𝑠𝑒 𝑛 ≠ 𝑚1, 𝑠𝑒 𝑛 = 𝑚
.
27
Com isso é possível definir a matriz métrica das funções de Hermite-Gauss,
([𝑮∅]𝑚𝑛), como,
⟨∅𝑚 , ∅𝑛⟩ = [𝑮∅]𝒎,𝒏 (3.1.3)
onde
𝑮∅ = [
1 0 ⋯ 00 1 ⋯ 0⋮ ⋮ ⋱ ⋮0 0 ⋯ 1
].
Ou seja, a métrica 𝑮∅ responde à matriz identidade.
3.2. Operador de derivação
A definição da matriz que representa o operador de derivação, assim como a matriz
métrica, varia com a função de base escolhida tendo em vista que cada função possui uma
derivada diferente. Sendo assim será necessário definir três operadores distintos, um
para cada uma das funções de base.
3.2.1. Operador de derivação de Chebyshev
Considere o operador �̂� =𝑑
𝑑𝑥. Se 𝐹 é o conjunto das funções contínuas em [-1,1], o
domínio de �̂� é um subconjunto (subespaço) de 𝐹, formado pelas funções diferenciáveis
em [-1,1], 𝐹1 = 𝐶1 ([−1,1]). Considere 𝑓 ∈ 𝐹1. Deve-se determinar a expansão de 𝑓 em
uma série envolvendo os polinômios de Chebyshev.
A expansão de 𝑓 é dada por
𝑓(𝑥) = ∑𝑓𝑘
∞
𝑘=0
𝑇𝑘(𝑥),
onde 𝑓𝑘 são constantes.
Aplicando o operador de derivação na expansão de 𝑓 obtém-se
(�̂�𝑓)(𝑥) =𝑑𝑓
𝑑𝑥= ∑𝑓𝑘
∞
𝑘=0
𝑑𝑇𝑘𝑑𝑥
(3.2.1.1)
28
Tomando o produto interno entre um polinômio de Chebyshev 𝑇𝑙 e (3.2.1.1),
obtém-se
⟨𝑇𝑙 , �̂�𝑓⟩ = ⟨𝑇𝑙,∑𝑓𝑘
∞
𝑘=0
𝑑𝑇𝑘𝑑𝑥⟩.
Como 𝑓𝑘 é constante,
=∑𝑓𝑘
∞
𝑘=0
⟨𝑇𝑙,𝑑𝑇𝑘𝑑𝑥
⟩. (3.2.1.2)
Define-se, então, a matriz �̃�𝒄 dada por:
⟨𝑇𝑙 ,𝑑𝑇𝑘𝑑𝑥
⟩ = [�̃�𝒄]𝑙𝑘
e, substituindo �̃�𝒄 em (3.2.1.1),
⟨𝑇𝑙 , �̂�𝑓⟩ = ∑(�̃�𝒄)𝑙𝑘𝑓𝑘
∞
𝑘=0
,
onde
∑(�̃�𝒄)𝑙𝑘𝑓𝑘
∞
𝑘=0
= (�̃�𝒄𝒇)𝑙.
Logo,
⟨𝑇𝑙, �̂�𝑓⟩ = (�̃�𝒄𝒇)𝑙
Os elementos da matriz �̃�𝒄 são definidos por
[�̃�𝒄]𝒍𝒌 = ⟨𝑇𝑙 ,𝑑𝑇𝑘𝑑𝑥
⟩ = {0, 𝑙 ≥ 𝑘
0, 𝑘 > 𝑙 𝑒 𝑘 + 𝑙 𝑝𝑎𝑟𝑘π, k > 𝑙 𝑒 𝑘 + 𝑙 ímpar
(3.2.1.3).
Com a definição (3.2.1.3) é possível construir a matriz para 0 ≤ 𝑙 ≤ 3 e 0 ≤ 𝑘 ≤ 3:
�̃�𝒄 = [
0 2π 0 4π0 0 3π 00 0 0 4π0 0 0 0
]
Por definição, obtêm-se a matriz 𝑫𝑪 que representa o operador de derivação dada
por,
29
�̃�𝒄 = 𝐆𝐂𝐃𝐂
onde que 𝐆𝐂 representa a matriz métrica.
3.2.2. Operador de derivação de Hermite
Para determinar a matriz que representa o operador de derivação dos polinômios
de Hermite, executa-se o mesmo processo descrito anteriormente com os polinômios de
Chebyshev.
Sabe-se que a derivada de um polinômio de Hermite de grau 𝑛,𝐻𝑛(𝑥), é dada por,
[6]
𝑑
𝑑𝑥𝐻𝑛(𝑥) = 2𝑛𝐻𝑛−1(𝑥) (3.2.2.1)
Tomando o produto interno entre a derivada de 𝐻𝑛(𝑥) e outro polinômio de
Hermite, 𝐻𝑚(𝑥), obtem-se:
[�̃�𝑯]𝒎,𝒏 =⟨𝐻𝑚 ,
𝑑
𝑑𝑥𝐻𝑛⟩ = ⟨𝐻𝑚 , 2𝑛𝐻𝑛−1⟩ = 2𝑛[𝐆
𝐇]𝒎,(𝒏−𝟏).
Assim,
[�̃�𝑯]𝒎𝒏 = 2𝑛[𝐆𝐇]𝒎,(𝒏−𝟏) (3.2.2.2).
Então �̃�𝑯 para 0 ≤ 𝑙 ≤ 2 e 0 ≤ 𝑘 ≤ 2 é dada por,
�̃�𝑯 = [0 11,31 00 0 67.880 0 0
],
Assim como para os polinômios de Chebyshev,
�̃�𝑯 = 𝑮𝑯𝑫𝑯
30
3.2.3. Operador de derivação de Hermite-Gauss
A derivada da função de Hermite-Gauss de grau n, ∅𝑛(𝑥), [7]
𝑑∅𝑛(𝑥)
𝑑𝑡= √
𝑛
2𝜎2∅𝑛−1(𝑥) − √
𝑛 + 1
2𝜎2∅𝑛+1(𝑥) (3.2.3.1)
Para definir a matriz que representa o operador de derivação, toma-se o seguinte
produto interno,
[�̃�∅]𝑙,𝑘 =⟨∅𝑙(𝑥),
𝑑∅𝑘(𝑥)
𝑑𝑥⟩ = √
𝑘
2𝜎2⟨∅𝑙(𝑥)∅𝑘−1(𝑥)⟩ − √
𝑘 + 1
2𝜎2⟨∅𝑙(𝑥)∅𝑘+1(𝑥)⟩
e tal produto interno permite definir
[�̃�∅]𝑙𝑘 =√𝑘
2𝜎2𝛿𝑙,𝑘−1 −√
𝑘 + 1
2𝜎2𝛿𝑙,𝑘+1 (3.2.3.2).
Dessa maneira, a matriz que representa o operador de derivação das funções de
Hermite-Gauss, possuem a diagonal principal nula enquanto as diagonais logo acima e
logo abaixo da principal, são diferentes de zero. Com a definição (3.2.3.2) é possível
construir a matriz �̃�∅ para 0 ≤ 𝑙 ≤ 2 e 0 ≤ 𝑘 ≤ 2:,
�̃�∅ =
[ 0 √
2
2𝜎20
√2
2𝜎20 √
3
2𝜎2
0 √3
2𝜎20
]
e
�̃�∅ = 𝑮∅𝑫∅
31
3.3. Operador Laplaciano (𝛁𝟐)
O operador Laplaciano aplicado a uma função unidimensional é
∇²𝑓 =𝑑2𝑓
𝑑𝑥2,
e esse operador pode também ser representado de forma matricial, do ponto de vista do
método espectral tau, assim como o operador de derivação. Essa representação será
apresentada a partir da matriz do operador de derivação e, dessa maneira, não será
necessário deduzi-la para cada função de base. Assim como o operador de derivação, cada
função de base possui sua matriz Laplaciano.
Para o caso dos polinômios de Chebyshev, o domínio do operador ∇² =𝑑2
𝑑𝑥2 é o
subconjunto (subespaço) das funções duas vezes diferenciáveis em [−1, 1], ℱ2 =
𝐶2([−1, 1]).
Considere 𝑓 ∈ ℱ2, 𝑔 ∈ ℱ e 𝑔(𝑥) =𝑑2𝑓
𝑑𝑥2 e que ℎ(𝑥) =
𝑑𝑓
𝑑𝑥, 𝑔(𝑥) =
𝑑ℎ
𝑑𝑥. As funções, ℎ(𝑥),
𝑔(𝑥) 𝑒 𝑓(𝑥), podem ser escritas por combinação linear de outras funções de base, 𝑏𝑘(𝑥),
de ordem k. Assim,
ℎ(𝑥) = ∑ℎ𝑘𝑏𝑘(𝑥)
𝑚
𝑘=0
; 𝑔(𝑥) = ∑𝑔𝑘𝑏𝑘(𝑥);
𝑚
𝑘=0
𝑓(𝑥) = ∑𝑓𝑘𝑏𝑘(𝑥)
𝑚
𝑘=0
Como foi descrito, ℎ(𝑥) é a derivada de 𝑓(𝑥), portanto,
ℎ(𝑥) =𝑑
𝑑𝑥𝑓(𝑥)
e
∑ℎ𝑘𝑏𝑘(𝑥)
𝑚
𝑘=0
=∑𝑓𝑘𝑑
𝑑𝑥𝑏𝑘(𝑥)
𝑚
𝑘=0
(3.3.1).
Tomando o produto interno da Equação (3.3.1) com outra função de base de grau
𝑙 obtém-se:
⟨𝑏𝑙(𝑥),∑ℎ𝑘𝑏𝑘(𝑥)
𝑚
𝑘=0
⟩ = ⟨𝑏𝑙(𝑥),∑𝑓𝑘𝑑
𝑑𝑥𝑏𝑘(𝑥)
𝑚
𝑘=0
⟩
∑ℎ𝑘⟨𝑏𝑙(𝑥), 𝑏𝑘(𝑥)⟩ = ∑𝑓𝑘 ⟨𝑏𝑙(𝑥),𝑑
𝑑𝑥𝑏𝑘(𝑥)⟩
𝑚
𝑘=0
𝑚
𝑘=0
(3.3.2)
De acordo com as definições de matriz métrica e de operadores de derivação, pode-
se inferir que a Equação (3.3.2) pode ser escrita como,
32
𝑮𝒉 = 𝒇𝑫 ̃ (3.3.3)
e, de (3.3.3), sabe-se que
𝒉 = 𝑮−𝟏�̃�𝒇 (3.3.4)
Utilizando o mesmo processo descrito mas agora para 𝑔(𝑥) =𝑑
𝑑𝑥ℎ(𝑥) obtém-se,
𝑮𝒈 = 𝑫 ̃𝒉 (3.3.5).
Substituindo a Equação (3.3.4) em (3.3.5) e isolando 𝒈𝒌, obtém-se,
𝒈 = 𝑮−𝟏𝑫 ̃𝑮−𝟏𝑫 ̃𝒇 = (𝑮−𝟏𝑫 ̃)𝟐𝒇 (3.3.6).
Define-se, então, a matriz 𝑫𝟐 como,
𝑫𝟐 = (𝑮−𝟏𝑫 ̃)𝟐
e, com isso, reescreve-se a Equação (3.3.6) como,
𝒈 = 𝑫𝟐𝒇 (3.3.7)
Está definido, então, pelas Equações (3.3.5) e (3.3.6) a matriz que representa o
operador Laplaciano, estabelecendo a Equação (3.3.8),
[𝑫𝟐]𝒍,𝒌 = ⟨𝑏𝑙,𝑑2
𝑑𝑥2𝑏𝑘⟩ (3.3.8)
33
Capítulo 4
Resultados e Validação do Método
Para validar o método espectral tau o problema escolhido foi a solução da equação
de Laplace em uma dimensão. A distribuição de potencial entre dois pontos, onde as
tensões estão determinadas, é descrita pela equação de Laplace e os potenciais nos pontos
são as condições de contorno. No Capítulo 4 desse trabalho são mostrados a solução
analítica da equação de Laplace unidimensional aplicada ao cálculo da distribuição de
potencial entre dois pontos, os resultados da implementação do método espectral tau e a
comparação entre a solução analítica e aproximada validando sua implementação. O
método espectral tau é implementado aqui utilizando as três funções de base estudadas,
os polinômios de Chebyshev, os polinômios de Hermite e as funções de Hermite-Gauss.
4.1. Equação de Laplace e o cálculo da distribuição de potencial entre dois pontos.
A solução analítica da distribuição do potencial entre dois pontos é simples. Esse
problema foi escolhido por se tratar de algo bastante conhecido que servirá como
validação da implementação do método espectral tau.
Considere os dois pontos com potenciais definidos na Figura (4-1).
Figura 4-1 - Geometria do problema. Pontos de potencial determinado.
34
g(x) é a função que representa a distribuição de potencial na direção x entre os pois
pontos. Os pontos estão separados por uma distância r e o potencial em x=-1 é
g(x = −1) = 0 e o potencial em x=1 é g(x = 1) = V.
O potencial entre os dois pontos se distribui respeitando a equação de Laplace,
onde ∇²g(x) = 0. Como o problema em questão é unidimensional, ou seja, varia ao longo
apenas de x, o operador ∇² se resume a
∇² =𝑑2
𝑑𝑥².
Dessa maneira,
∇2g(x) =𝑑²g
𝑑𝑥² (4.1.1)
e a solução da equação diferencial anterior determina que a função g(x) é dada por [11]
g(x) =V
rx (4.1.2)
Supondo que V=2 e r=2, posicionando as condições de contorno em 𝑥 = −1 e em
𝑥 = 1, obtem-se que
g(x) = x + 1
e a distribuição do potencial pode ser representado pela a figura a seguir:
Figura 4-2 - Distribuição de potencial entre dois pontos
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 10
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
Posição dos pontos (x)
Pote
ncia
l (V
)
Distribuição do potencial entre dois pontos
35
4.2. Método espectral tau aplicado à distribuição de potencial entre dois pontos
Para a implementação do método espectral tau foi desenvolvida uma ferramenta
computacional baseada no software MATLAB® por se tratar de um excelente programa
para problemas que envolvam matrizes. Algumas funções, como a função que calcula o
valor do polinômio de Chebyshev no ponto especificado, não foram desenvolvidas nesse
trabalho e foram obtidas em pesquisas externas. As funções desenvolvidas são
responsáveis por calcular as matrizes métricas 𝑮, as matrizes de derivada 𝑫, as matrizes
do Laplaciano 𝑫² e os valores das funções de base nos determinados pontos.
Para aplicar o método espectral tau é necessário definir a modelagem matemática
corretamente para sua implementação. A geometria do problema será limitada pelos
polinômios de Chebyshev tendo em vista que tais funções são operáveis apenas entre -1
e +1. Definindo as condições de contorno nas fronteiras do domínio dos polinômios de
Chebyshev, pode-se definir as condições de contorno do problema como
g(x = −1) = 0
g(x = 1) = 20.
Como definido anteriormente, deve-se adicionar modos para as condições de
contorno. Como temos duas condições de contorno, serão necessários dois modos à mais
que serão definidos pela variável "𝑁𝐶𝐶”.
A equação diferencial submetida ao método espectral tau é a equação de Laplace
definida pela Equação (4.1.1). Então deve-se aproximar a função ϕ(x) por uma
combinação linear de uma função de base ortogonal. Assim
ϕ(x) ≅ f(x) = ∑𝑢𝑘𝑏𝑘
𝑛
𝑘=0
(𝑥)
onde 𝑛 representa os modos necessários para acertar a parte diferencial do
sistema linear proposto. O número total de modos é definido por
𝑚 = 𝑛 + 𝑁𝐶𝐶 (4.2.1)
O somatório que antes ia de 0 à 𝑛 deve agora ir até 𝑚 e deve ser diferenciado duas
vezes. Assim,
𝑑2𝑓(𝑥)
𝑑𝑥²=𝑑2
𝑑𝑥2∑𝑢𝑘𝑏𝑘(x) = 0
𝑚
𝑘=0
(4.2.2)
36
De acordo com as condições de contorno impostas, a função aproximada deve
também respeitá-las. Assim,
{
𝑑
2𝑓(𝑥)
𝑑𝑥²=𝑑2
𝑑𝑥2∑𝑢𝑘𝑏𝑘(x) = 0
𝑚
𝑘=0
𝑓(−1) = ∑𝑢𝑘𝑏𝑘(−1) = 0
𝑚
𝑘=0
𝑓(1) = ∑𝑢𝑘𝑏𝑘(1) = 𝑉
𝑚
𝑘=0
(4.2.3)
Tomando o produto interno entre a primeira equação do Sistema (4.2.3) e outra
função de base de ordem "𝑙" onde 0 ≤ 𝑙 ≤ 𝑚 obtém-se
{
⟨𝑏𝑙(𝑥),
𝑑2
𝑑𝑥2∑𝑢𝑘𝑏𝑘(x)
𝑚
𝑘=0
⟩ = 0
∑𝑢𝑘𝑏𝑘(−1)
𝑚
𝑘=0
= 0
∑𝑢𝑘𝑏𝑘(1)
𝑚
𝑘=0
= 𝑉
(4.2.4)
Como as incógnitas 𝑢𝑘 são constantes, elas saem do operador Laplaciano e
também do produto interno. Assim o Sistema (4.2.4) torna-se
{
∑𝑢𝑘
𝑚
𝑘=0
⟨𝑏𝑙(𝑥),𝑑2
𝑑𝑥2bk(𝑥)⟩ = 0
∑𝑢𝑘𝑏𝑘(−1)
𝑚
𝑘=0
= 0
∑𝑢𝑘𝑏𝑘(1)
𝑚
𝑘=0
= 𝑉
(4.2.5)
A primeira equação do Sistema (4.2.5) pode ser reescrita de maneira matricial da
seguinte maneira,
∑𝑢𝑘
𝑚
𝑘=0
⟨𝑏𝑙(𝑥),𝑑2
𝑑𝑥2bk(𝑥)⟩ = 0
[ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2b0(𝑥)⟩ ⋯ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2bm−1(𝑥)⟩ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2bm(𝑥)⟩
⋮ ⋱ ⋮ ⋮
⟨𝑏𝑚−1(𝑥),𝑑2
𝑑𝑥2b0(𝑥)⟩ ⋯ ⟨𝑏𝑚−1(𝑥),
𝑑2
𝑑𝑥2bm−1(𝑥)⟩ ⟨𝑏𝑚−1(𝑥),
𝑑2
𝑑𝑥2bm(𝑥)⟩
⟨𝑏𝑚(𝑥),𝑑2
𝑑𝑥2b0(𝑥)⟩ ⋯ ⟨𝑏𝑚(𝑥),
𝑑2
𝑑𝑥2bm−1(𝑥)⟩ ⟨𝑏𝑚(𝑥),
𝑑2
𝑑𝑥2bm(𝑥)⟩ ]
[
𝑢0⋮
𝑢𝑚−1𝑢𝑚
] = [
0⋮00
]
37
e que, de maneira geral, pode ser escrita como
[𝑫𝟐]𝑚,𝑚[𝑼]𝑚,1 = [𝟎]𝑚,1 (4.2.6)
onde [𝑼]𝑚,1 é uma matriz coluna com 𝑚 linhas onde estão as constantes 𝑢𝑘, 𝑫𝟐 é a matriz
que representa o operador Laplaciano definida pela Equação (3.3.8) e [𝟎]𝒎,𝟏 uma matriz
coluna de zeros, que assim como [𝑼]𝑚,1, possui 𝑚 linhas.
Mas, de acordo com a definição do método espectral tau, 𝑚 = 𝑛 + 𝑁𝐶𝐶, onde 𝑁𝐶𝐶
são os modos adicionais incluídos à aproximação, para respeitar as condições de
contorno. Como o problema envolve duas condições de contorno, 𝑁𝐶𝐶 = 2 e, para ajustar
o Sistema (4.2.6), é necessário, então, realizar a substituição das suas duas últimas linhas
pelas equações das condições de contorno. O Sistema (4.2.7) reformulado é
[ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2b0(𝑥)⟩ ⋯ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2bm−1(𝑥)⟩ ⟨𝑏0(𝑥),
𝑑2
𝑑𝑥2bm(𝑥)⟩
⋮ ⋱ ⋮ ⋮b0(−1) ⋯ bm−1(−1) bm(−1)
b0(1) ⋯ bm−1(1) bm(1) ]
[
𝑢0⋮
𝑢𝑚−1𝑢𝑚
] = [
0⋮0𝑉
] (4.2.7)
e que pode ser reescrito como
𝑨𝑼 = 𝒃,
onde a matriz 𝑨 é a antiga matriz de Laplace do Sistema (4.2.6), 𝑫𝟐, com as duas últimas
linhas substituídas, assim como o vetor 𝒃 é o antigo vetor de zeros , 𝟎, também com as
duas últimas linhas substituídas. A matriz 𝑨 pode ser facilmente determinada tendo em
vista que de acordo com a Equação (3.4.1) pode-se obter a matriz Laplaciano e, de acordo
com as Equações (3.1.1) (3.1.2) e (3.1.3), é possível obter os produtos internos entre as
funções de base.
Suprimindo os índices do Sistema (4.2.7) e multiplicando os dois lados pela inversa
de 𝑨, obtém-se
𝑼 = 𝑨−𝟏𝒃 (4.2.8)
No caso estudado, a matriz 𝑨 apresentou mal condicionamento identificado pelo
software. A solução dada foi implementar uma aproximação por mínimos quadrados
onde,
𝑨−𝟏 ≅ (𝑨𝑻𝑨)−𝟏𝑨𝑻.
e, obtendo os coeficientes 𝑢𝑘, é possível construir a função 𝑓(𝑥) uma vez que
𝑓(x) = ∑𝑢𝑘𝑏𝑘
𝑚
𝑘=0
(𝑥).
38
De posse da função 𝑓(𝑥) obtida para cada uma das três funções de base, os
polinômios de Chebyshev, os polinômios de Hermite e as funçõs de Hermite-Gauss, é
possível compará-la com a solução analítica do problema discutida no Tópico 4.1.
4.2.1. Resultados obtidos com os polinômios de Chebyshev
Escolhendo o valor de 𝑚 como 10, ou seja, com 10 modos totais, sendo 8 para a
parte diferencial do Sistema (4.2.6) e 2 modos para atender as condições de contorno, foi
possível implementar o método espectral tau. Plotando as soluções analítica e
aproximada em mesmo gráfico, obteve-se a Figura 4-3 a seguir
Figura 4-3 - Distribuição de potencial entre dois pontos (polinômios de Chebyshev)
4.2.2. Resultados obtidos com os polinômios de Hermite
Assim como para os polinômios de Chebyshev, escolhendo o valor de 𝑚 como 10,
ou seja, com 10 modos totais, sendo 8 para a parte diferencial do Sistema (4.2.6) e 2 modos
para atender as condições de contorno, foi possível implementar o método espectral tau.
Plotando as soluções analítica e aproximada em mesmo gráfico, obteve-se a Figura 4-4 a
seguir
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 10
2
4
6
8
10
12
14
16
18
20
Posição dos pontos (x)
Pote
ncia
l (V
)
Distribuição de potencial entre dois pontos em -1 e 1. Analitica e Aproximada
Solução Analitica
Solução aproximada
39
Figura 4-4 - Distribuição de potencial entre dois pontos (polinômios de Hermite)
4.2.3. Resultados obtidos com as funções de Hermite-Gauss
Assim como para os polinômios de Chebyshev e Hermite, escolhendo o valor de 𝝈
igual a 1, 𝑚 como 10, ou seja, com 10 modos totais, sendo 8 para a parte diferencial do
Sistema (4.2.6) e 2 modos para atender as condições de contorno, foi possível
implementar o método espectral tau. Plotando as soluções analítica e aproximada em
mesmo gráfico, obteve-se a Figura 4-5 a seguir
Figura 4-5 - Distribuição de potencial entre dois pontos, 𝜎 = 1 (funções de Hermite-Gauss)
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 10
2
4
6
8
10
12
14
16
18
20
Posição dos pontos (x)
Pote
ncia
l (V
)
Distribuição de potencial entre dois pontos (Hermite)
Solução Aproximada
Solução Analitica
-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1-150
-100
-50
0
50
100
Posição dos pontos (x)
Pote
ncia
l(V
)
Distribuição de potencial entre dois pontos (Hermite-Gauss)
Solução Aproximada
Solução Analitica
40
Aumentando o valor de 𝝈 para 10 e mantendo os outros parâmetros iguais, obteve-
se a Figura 4-6
Figura 4-6 - Distribuição de potencial entre dois pontos, 𝜎 = 10 (funções de Hermite-Gauss)
4.3. Validação do método
Observando a Figura (4-3) e Figura (4-4), pode-se dizer que o método espectral
tau, como previu Hermann [10], aproximou com sucesso a solução da Equação de Laplace
pois, de maneira geral, as soluções aproximadas seguiram a forma da solução analítica.
Nestes casos, as funções de base escolhidas foram os polinômios de Chebyshev e os
polinômios de Hermite.
Na implementação do método para as funções de Hermite-Gauss houve, a
princípio, uma discrepância grande entre os valores analíticos e aproximados de
potencial. Essa diferença ocorreu para o menor valor da largura da Gaussiana (𝜎). Para 𝜎
igual a 10, a função convergiu para o resultado analítico e quanto maior for o valor de 𝜎,
menor a diferença entre a solução analítica e a aproximada. Isso acontece porque quanto
maior for o valor de 𝜎, maior será a largura da Gaussiana fazendo com que sua
convergência para 0 se afaste do domínio do problema (entre -1 e 1). Desta maneira, no
domínio especificado, as funções de Hemite-Gauss tendem a se comportar como os
próprios polinômios de Hermite para valores elevados de 𝜎, que como visto na Figura 4-
4, aproximaram a solução analítica com sucesso.
41
Capítulo 5
Conclusão de Perspectivas
O método espectral tau aplicado à equação de Laplace Unidimensional serviu para
mostrar que o método é realmente válido para resolver problemas com equações
diferenciais parciais e problemas de valor de contorno. A sua validação ocorreu a partir
do momento que não houve grande diferença entre os valores obtidos na aproximação da
distribuição de potencial e a sua solução analítica.
Devido à correta caracterização do problema da eletrostática escolhido e da
correta modelagem matricial do método espectral tau e das funções de base, as
implementações numéricas foram obtidas com rapidez no ambiente MATLAB®.
O programa computacional que implementa o método no software MATLAB® foi
construído de maneira que as funções, construídas separadamente, podem ser alteradas
com facilidade. Assim o programa ficou limpo e de fácil compreensão.
Para trabalhos futuros pretende-se aplicar o método espectral tau em outros
problemas de valor de contorno e problemas que envolvam evolução temporal como a
equação de onda. Coloca-se também em pauta, a procura de uma explicação para o mal
condicionamento da matriz 𝑨.
42
Referências Bibliográficas
[1] ORSZAG, GOTTLIEB "Numerical Analysis of Spectral Methods" Theory and
Applications Book.
[2] JÁUBER O. C. "Métodos espectrais e seus derivados", Universidade Federal
de Santa Catarina, Florianópolis, SC, Brasil
[3] GARDNER D. R., TROGDON S. A., DOUGLASS R. W. "A modified Tau Spectral
Method That Eliminates Spurius Eingenvalues"
[4] MASON J. C., HANDSCOMB D.C. "Chebyshev Polynomials" book 2003
[5] WEISSTEIN, ERIC W. "Chebyshev Polynomial of the First Kind."
From MathWorld--A Wolfram Web Resource.
http://mathworld.wolfram.com/ChebyshevPolynomialoftheFirstKind.html
Data de acesso: 10/08/2014
[6] FURTADO M. F., "Algumas Realizações de Charles Hermite" Universidade de
Brasilia, 1996.
[7] BAYIN S. S., "Mathematical Methods in Science and Engineering" New York,
2006.
[8] LANCZOS C. and ORTIZ E. L., "Step by step tau method - part I. Piecewise
polynomial approximations", Imperial College, University of London, London, England.
[9] FOX L. and PARKER I. B., "Chebyshev polynomials in numerical analysis"
Oxford University Press, Oxford, UK 1968.
[10]HERMANN R. "Spectral Methods for Partial Differential Equations"
Engineering Sciences and Applied Mathematics.
[11] FONTANA E. "Eletromanetismo - Parte 1 - Problemas de Valores de
Fronteira em Eletrostática" Edição 01.2011, Departamento de Eletrônica e Sistemas
UFPE.
[12] ANANTHASAYANAM B. "ortho_poly.m - Compute orthogonal
polynomials" 2005. http://ceta.mit.edu/comp_spec_func/ Data de acesso: 12/08/14.
43
[13] SUINESIAPUTRA A. "hermite.m - compute the hermite polynomials"
2010http://www.mathworks.com/matlabcentral/fileexchange/27746-hermite-
polynomials/content/hermite.m Data de acesso: 12/08/14
44
Anexos
Códigos MATLAB
Em anexo seguem os códigos comentados utilizados na implementação do método
espectal tau. A explicação sobre o código da implementação em si, está descrita com os
polinômios de Chebyshev. Para implementar as outras funções de base, basta trocar as
funções que calculam o valor nos pontos (ortho_poly, Hermite, Hgauss3), trocar a matriz
métrica (métrica, metricaH, e metricaHG) e as matrizes de derivação dentro de “lnp”
(derivada, derivadaH, derivaHG).
A função ortho_poly.m que calcula o valor dos polinômios de Chebyshev não foi
elaborada pelo o autor desse trabalho e está disponível para livre distribuição[12], assim
como a função hermite.m[13].
A-1 Laplace 1D
Este script calcula numericamente a solução da equação de Laplace em uma
dimensão usando o método espectral tau. Calcula-se a distribuição de potencial entre dois
pontos localizados em x = -1 e x = 1.
O número total de modos e mx+1. Os primeiros nx+1 modos são usados ajustar a
função incógnita a partir da EDP. Um sistema de nx+1 equações lineares com mx+1
variáveis são obtidas. Os demais ncc = mx-nx modos são usados para ajustar as
condições de contorno. NCC é o número de pontos escolhidos sobre a fronteira para
ajustar as condições de contorno.
clear all
clc
mx = 2;
nx = 0;
ncc = mx - nx;
45
Determina a matriz A que está na parte direita do sistema de equações lineares
Au = 0. Esta matriz é quadrada de ordem (mx+1). Apenas (nx+1) linhas desta matriz
serão aproveitadas. As demais devem ser ajustadas para incorporar as condições de
contorno.
A = lnp(mx,0);
O vetor x contém as coordenadas da fronteira sobre as quais as condições de
contorno são medidas.
x = linspace(-1,1,ncc);
As linhas da matriz A que correspondem as equações algébricas referentes as
condições de contorno devem ser ajustadas para satisfazer as CC. No loop a seguir, este
ajuste e efetuado.
b = zeros(mx+1,1);
cc = 1;
for i1 = 0:mx
k = i1 + 1; %k identifica a linha da matriz A a ser ajustada.
if (i1 > nx)
for i2 = 0:mx %l identifica a coluna da matriz A.
l = i2 + 1; %Calcula o valor da função de base de grau i2 no ponto
%x(cc), ou seja, em determinado ponto no contorno.
A(k,l) = ortho_poly(1,x(cc),i2);
end
cc = cc + 1;
end
end
Agora, determina os valores de contorno em cada par coordenado sobre a
fronteira.
cc = 1;
for i1 = 0:mx
k = i1 + 1; %k identifica a linha da matriz b.
if (i1 > nx)
b(k, 1) = ft(x(cc),0,0);
46
cc = cc + 1;
end
end
u = inv(A)*b;
X = linspace(-1,1,10);
pot = zeros(1,length(X));
for i1 = 0:mx
k = i1 + 1;
%k identifica a linha da matriz u.
pot = pot + u(k)*ortho_poly(1,X,i1);
end
plot(X,pot,'*k')
hold on
X1 = linspace(-1,1,15);
A função fi eh a solução analítica da equacao de Laplace para avaliar a distribuição
de potencial entre dois pontos onde a ddp eh 20.
fi=10.*X1+10;
plot(X1,fi,'+r')
A-2 Lnp.m
A função lnp.m é usada para calcular a matriz que representa o operador
Laplaciano.
function[L_]=lnp(nx,ny) %calcula a matriz que representa o operador
Laplaciano
Dx_=derivada(nx); %a derivada,m pode ser mudada para qualquer função de
base.
L_=Dx_^2;
End
47
A-3 Derivada.m
A função derivada.m é usada para calcular a matriz que representa o operador de
derivação dos polinômios de Chebyshev.
function[D_]=derivada(n)
D_=zeros(n+1,n+1);
for l=1:(n+1);
for k=(l+1):(n+1);
q=mod(k+l,2);%q=1 se k+l for impar, q=0 se k+l for par
if q==1
D_(l,k)=k*pi;
end
end
end
end
A-4 DerivadaH.m
A função derivadaH.m é usada para calcular a matriz que representa o operador de
derivação dos polinômios de Hermite.
for i=1:n+1;
for j=1:n+1;
if i==(j-1)
D_(i,j)=sqrt(2)*2*j*(2^(j-1))*factorial(j-1);
end
end
end
end
48
A-5 DerivaHG
A função derivaHG.m é usada para calcular a matriz que representa o operador de
derivação das funções de Hermite-Gauss.
function [D]=derivaHG(sigma,l,k)
for i=1:l+1;
for j=1:k+1;
p1=sqrt(j/(2*sigma^2));
p2=sqrt((j+1)/(2*sigma^2));
if i==j-1
p11(i,j)=p1;
else
p11(i,j)=0;
end
if i==j+1
p21(i,j)=p2;
else
p21(i,j)=0;
end
D(i,j)=p11(i,j)+p21(i,j);
end
end
end
A-6 Metrica.m
A função métrica.m calcula a matriz métrica que representa o produto interno
entre polinômios de Chebyshev.
function[G]=metrica(n,m)
G = zeros(n+1,m+1);
for i=1:(n+1);
for j=1 : (m+1);
49
if i==j
G(i,j)=sqrt(pi)*(2^i)*factorial(i);
end
end
end
end
A-7 metricaH.m
A função métricaH.m calcula a matriz métrica que representa o produto interno
entre polinômios de Hermite.
G = zeros(n+1,m+1);
for i=1:(n+1);
for j=1 : (m+1);
if i==j
G(i,j)=sqrt(pi)*(2^i)*factorial(i);
end
end
end
end
A-8 metricaHG.m
A função métricaHG.m calcula a matriz métrica que representa o produto interno
entre funções de Hermite-Gauss.
function [M]=metricaHG(l,k)
for i=1:l+1;
for j=1:k+1;
if i==j
M(i,j)=1;
else
M(i,j)=0;
50
end
end
end
A-9 ortho_poly.m
A função ortho_poly.m calcula o valor do polinômio de Chebyshev do tipo “kf”, de
grau “n” nos valores de “x”. Essa função não foi elaborada pelo autor desse trabalho.
function pl=ortho_poly(kf,x,n)
% This is a code downloaded from the website of MIT.
% http://ceta.mit.edu/comp_spec_func/
%
==========================================================
% Purpose: Compute orthogonal polynomials: Tn(x) or Un(x),
% or Ln(x) or Hn(x), and their derivatives
% Input : KF --- Function code
% KF=1 for Chebyshev polynomial (First kind) Tn(x)
% KF=2 for Chebyshev polynomial (Second kind) Un(x)
% n --- Order of orthogonal polynomials
% x --- Argument of orthogonal polynomials
% Output: PL(n) --- Tn(x) or Un(x) or Ln(x) or Hn(x)
% DPL(n)--- Tn'(x) or Un'(x) or Ln'(x) or Hn'(x)
%
=========================================================
% The only improvement in this program is it accepts vector arguments for
x
% make sure that x is a row or column vector and not a matrix.
[r,c]=size(x);
if r==1 | c==1
rowvec = 0;
51
if r==1
x=x';
rowvec = 1;
end
else
error('x must be a vector, and cannot be a matrix');
end
lenx = length(x);
if n==0
if rowvec
pl = ones(1,lenx);
else
pl = ones(lenx,1);
end
else
pl = zeros(lenx,n);
a=2;
b=0;
c=1;
y0=1;
y1=2.*x;
% the i'th position in pl corresponds to the i'th term
% don't bother storing pl = 1;
pl(:,1)=2.*x;
if (kf == 1)
y1=x;
pl(:,1)=y1;
end
52
for k=2:n
yn=(a.*x+b).*y1-c*y0;
pl(:,k)=yn;
y0=y1;
y1=yn;
end
if rowvec
pl = pl(:,n)';
else
pl = pl(:,n);
end
end
A-10 Hermite.m
A função Hermite.m calcula o valor do polinômio de Hermite de grau “n” nos
valores de “x”. Essa função não foi elaborada pelo autor desse trabalho.
function h = hermite(n,x)
% HERMITE: compute the Hermite polynomials.
%
% h = hermite(n)
% h = hermite(n,x)
%
% Inputs:
% - n is the order of the Hermite polynomial (n>=0).
% - x is (optional) values to be evaluated on the resulting Hermite
% polynomial function.
%
% There are two possible outputs:
% 1. If x is omitted then h is an array with (n+1) elements that contains
% coefficients of each Hermite polynomial term.
53
% E.g. calling h = hermite(3)
% will result h = [8 0 -12 0], i.e. 8x^3 - 12x
%
% 2. If x is given, then h = Hn(x) and h is the same size of x.
% E.g., H2(x) = 4x^2 - 2
% calling h = hermite(2,[0 1 2])
% will result h = [-2 2 14]
%
% More information:
% - about the Hermite polynomial:
http://mathworld.wolfram.com/HermitePolynomial.html
% - some examples of this function:
% http://suinotes.wordpress.com/2010/05/26/hermite-polynomials-with-
matlab/
%
% Authors: Avan Suinesiaputra ([email protected])
% Fadillah Z Tala ([email protected])
% rev.
% 26/05/2010 - first creation.
% - bug fixed: error when hermite(0,x) is called (x isn't empty)
% 24/09/2010 - bug fixed: the size of x does match with y in line 50.
% (thanks to Shiguo Peng)
% check n
if( n<0 ), error('The order of Hermite polynomial must be greater than or
equal to 0.'); end
% again check n is an integer
if( 0~=n-fix(n) ), error('The order of Hermite polynomial must be an
integer.'); end
% call the hermite recursive function.
54
h = hermite_rec(n);
% evaluate the hermite polynomial function, given x
if( nargin==2 )
y = h(end) * ones(size(x));
p = 1;
for i=length(h)-1:-1:1
y = y + h(i) * x.^p;
p = p+1;
end
% restore the shape of y, the same as x
h = reshape(y,size(x));
end
function h = hermite_rec(n)
% This is the reccurence construction of a Hermite polynomial, i.e.:
% H0(x) = 1
% H1(x) = 2x
% H[n+1](x) = 2x Hn(x) - 2n H[n-1](x)
if( 0==n ), h = 1;
elseif( 1==n ), h = [2 0];
else
h1 = zeros(1,n+1);
h1(1:n) = 2*hermite_rec(n-1);
h2 = zeros(1,n+1);
h2(3:end) = 2*(n-1)*hermite_rec(n-2);
h = h1 - h2;
55
end
A-11 Hgauss3.m
A função Hgauss3.m calcula o valor da função de Hermite-Gauss com largura
“sigma”, de ordem “n”, nos valores de t.
function [Y]=Hgauss3(sigma,t,n)
x1=t./sigma;
parte1=1/(sqrt(2^n*factorial(n)*sigma*sqrt(pi)));
parte2=exp(-1.*((x1.^2)./2));
Y=parte1.*parte2.*hermite(n,x1);
end
56
