Expansão em caos polinomial

Introducao e toolbox ECP em EDE Estudo de caso: EDO estocastica Estudo de caso: Oscilador aleatorio Referencias

Expansao em caos polinomial

Wilson N. de Freitas

Departamento de Engenharia Eletrica — PUC–Rio

31 de Agosto de 2007


1 Introducao e toolboxExpansao em Caos Polinomial (ECP)Espaco de HilbertExpansoes ortogonaisEspaco das funcoes quadraticamente integraveis L2(D)Polinomios ortogonaisDefinicao de Caos Polinomial

2 ECP em EDEEquacoes diferenciais estocasticasMetodo de Galerkin

3 Estudo de caso: EDO estocasticaEDO com termo aleatorio

4 Estudo de caso: Oscilador aleatorioOscilador aleatorio de segunda ordem

5 Referencias


Expansao em Caos Polinomial (ECP)

Caos Polinomial

Proposo por Norbert Wiener em 1938

Emprega polinomios de Hermite em termos de variaveisaleatorias Gaussianas para descrever variaveis aleatorias.

E como e que isso acontece?


Espaco de Hilbert

Espaco de Hilbert H

e um espaco vetorial sobre um corpo F (que pode ser R ou C)

possui produto interno 〈·, ·〉

e completo como espaco metrico, com relacao a metrica (d(·, ·))gerada pela norma (‖ · ‖) induzida pelo produto interno

‖v‖ =√〈v, v〉, onde v ∈ H

ed(u, v) = ‖u− v‖, onde u, v ∈ H

Espacos de Hilbert populares

(Rn; 〈·, ·〉)

(Cn; 〈·, ·〉)

(L2(D); 〈·, ·〉)


Expansoes ortogonais

Definition

Um conjunto de vetores Φ ∈ H e um conjunto ortonormal se qualquerpar de vetores distintos φi, φj ∈ Φ sao ortogonais, isto e, 〈φi, φj〉 = 0sempre que i 6= j e adicionalmente, ‖φi‖ = 1 para cada φi ∈ Φ

Theorem

Teorema das series de Fourier: Seja Φ = φnn∈N um conjuntoortonormal contavel em um espaco de Hilbert H, entao, Φ e uma baseortonormal de H. Cada y ∈ H tem uma unica expansao em Φ

y =∑n∈N

〈y, φn〉φn


Espaco das funcoes quadraticamente integraveis L2(D)

L2(D)

L2(D) e o espaco das funcoes quadraticamente integraveis em umdomınio D

f ∈ L2(D) se

∫D

|f(x)|2dx < ∞

O produto interno em L2(D) e

〈f, g〉 =

∫D

f(x)g(x)dx

Variaveis aleatorias com variancia finita tambem pertencem a L2(D)

E[|X|2

]=

∫D

|x|2dP (x) =

∫D

|x|2f(x)dx < ∞

Nesse caso o produto interno e

〈X, Y 〉 = E[XY

]=

∫D

xydP (x, y) =

∫D

xyf(x, y)dxdy


Polinomios ortogonais


Seja Qn(x)∞n=0 um conjunto de polinomios e n e o grau do polinomio

Quando os polinomios sao funcoes de uma variavel aleatoria X temos aseguinte relacao

E[Qn(X)Qm(X)

]=

∫D

Qn(x)Qm(x)dP (x)

=

∫D

Qn(x)Qm(x)f(x)dx

= hnδnm

Variaveis aleatorias Polinomios Domınio

Gaussiana Hermite (−∞,∞)Gama Laguerre [0,∞)Beta Jacobi [a, b]

Uniform Legendre [a, b]Poisson Charlier 0, 1, . . . Binomial Krawtchouk 0, 1, . . . , N




Polinomios de Hermite Hn(x)

Definicao e−x2Hn(x) = (−1)n dn

dxn

(e−x2)

Ortogonalidade 1√π

∫∞−∞ e−x2

Hm(x)Hn(x)dx = 2nn!δmn

Polinomios de Laguerre L(α)n (x)

Definicao e−xxαL(α)n (x) = 1

n!dn

dxn

(e−xxn+α

), α > −1

Ortogonalidade∫∞0 e−xxαL

(α)m (x)L

(α)n (x)dx =

Γ(n+α+1)n! δmn

Polinomios de Charlier Cn(x; a)

Definicao ax

x! Cn(x; a) = ∇n(

ax

x!

), a > 0

Ortogonalidade∑∞

x=0ax

x! Cm(x; a)Cn(x; a) = a−nean!δmn


Definicao de Caos Polinomial


Considere Θ o conjunto de todas as variaveis aleatorias com varianciafinita relacionadas ao espaco amostral Ω. Se ξ ∈ Θ entao ξ : Ω → R.

Para cada ξ ∈ Θ entao∫

D|ξ|2dP (ξ) < ∞ e portanto Θ e um espaco de

Hilbert.

Seja Φp∞p=0 um conjunto de polinomios ortogonais em Θ, logo,Φp : Θ → Θ.

Pelo teorema de series de Fourier pode-se expandir qualquer elemento deΘ em Φp∞p=0

X(ω) =

∞∑i=0

xiΦi(ξ(ω))

O conjunto Φp∞p=0 e o Caos Polinomial



O que foi visto ate agora?

Definicao de espacos de Hilbert

Teorema de series de Fourier

L2(D)Polinomios ortogonais

E o que fazer com tudo isso ?


Equacoes diferenciais estocasticas


Considere a equacao diferencial estocastica:

Λu = f

Λ ≡ Λ(x, t, ω) e um operador diferencial estocastico comderivadas em t e x e com um termo aleatorio ω

u ≡ u(x, t;ω) e a funcao incognita

x e t sao as variaveis independentes

f ≡ f(x, t;ω) e o termo de excitacao, aleatoria ou nao




A solucao u pode ser expandida em caos polinomial

u(x, t; ω) =

∞∑i=0

ui(x, t)Φi(ξ(ω))

Para obter uma aproximacao analıtica da solucao e necessario truncar aserie em um numero P finito de termos.

u(x, t; ω) =

P∑i=0

ui(x, t)Φi(ξ(ω))

O truncamento introduz um erro de aproximacao na solucao u.

Considere o erro r(x, t) de aproximacao como

r(x, t) = Λu− f

quando u e representado pela serie truncada.


Metodo de Galerkin

Metodo de Galerkin

Exigindo que o erro r seja ortogonal ao subespaco gerado peloconjunto ΦiP

i=0 temos que:

〈r(x, t),Φi〉 = 0

onde i = 0, 1, . . . , P

Com isso〈Λu− f,Φi〉 = 0

〈Λu, Φi〉 = 〈f,Φi〉〈Λ(

∑j ujΦj),Φi〉 = 〈f,Φi〉∑

j Λuj〈Φj ,Φi〉 = 〈f,Φi〉Λui = 〈f,Φi〉

temos um conjunto de P equacoes diferenciais deterministicas.

Esta abordagem e conhecida como metodo de Galerkin.


EDO com termo aleatorio


Consideremos a seguinte equacao diferencial ordinaria com termoaleatorio

dy(t)dt

= −ky(t)

k ≡ k(ω) = k(ξ(ω)) e uma variavel aleatoria com funcao deprobabilidade f(k)

A solucao y(t) pode ser redefinida como y(t, ω) = y(t, ξ(ω)) eportanto pode ser aplicada a ECP.

y(t, ω) =P∑

i=0

yi(t)Φi(ξ(ω))

A variavel aleatoria k(ω) pode ser escrita na mesma base de y(t, ω)

k(ω) =P∑

i=0

kiΦi(ξ(ω))



Aplicando a ECP

Substituindo as expansoes na equacao diferencial temos:

P∑i=0

dy(t)dt

Φi = −P∑

i=0

P∑j=0

ΦiΦjkiyj(t)

Aplicando o metodo de Galerkin a equacao acima temos:

dyl(t)dt

〈Φl,Φl〉 = −P∑

i=0

P∑j=0

〈ΦiΦj ,Φl〉kiyj(t)

onde l = 0, 1, . . . , P .

Este sistema de equacoes diferenciais deterministicas pode serresolvido com metodos numerios convencionais, como exemplo:metodo de Euler ou Runge-Kutta de segunda ou quarta ordem.Ainda faltam as condicoes iniciais.



Valor esperado da solucao

Tomando o valor esperado da solucao y(t, ω) na ECP

E[y(t, ω)

]=

P∑i=0

yi(t)E[Φi(ξ(ω))

]O valor esperado esta intimamente ligada ao produto interno nabase Φi∞i=0, logo:

E[Φi(ξ(ω))

]= 〈Φi(ξ(ω)),Φ0(ξ(ω))〉 = 0

para todo i > 0, pois Φ0 = 1 para a maioria dos casos.

O valor esperado de uma ECP sempre recai sobre o seu primeirotermo.

E[y(t, ω)

]= y0(t)



Condicoes de contorno

As condicoes de contorno recaem sobre o primeiro termo da ECP,como foi visto no metodo de Galerkin.

E[y(0, ω)

]= y0(0) = y0 = constante

De posse das condicoes iniciais e possıvel resolver o sistema deequacoes diferenciais

dyl(t)dt

〈Φl,Φl〉 = −P∑

i=0

P∑j=0

〈ΦiΦj ,Φl〉kiyj(t)

para l = 0, 1, . . . , P .

O objetivo e encontrar a y0(t) numericamente, pois este termorepresenta o valor esperado da solucao y(t, ω).



Resolvendo a EDO sem a ECP

A condicao inicialy(0, ω) = y0

constante em qualquer cenario. Tambem poderia ser uma variavelaleatoria, contanto que fosse ortogonal a k para nao complicar nas contas.

O valor esperado da solucao estocastica e dado por:

E[y(t, ω)

]= y0

∫D

e−ktf(k)dk

De posse desta solucao analıtica pode-se comparar com a ECP.

Considerando o seguinte erro

ε(t) =∣∣∣ 〈yECP (t)〉 − 〈y(t)〉

〈y(t)〉

∣∣∣onde 〈y(t)〉 e o valor esperado da solucao estocastica e〈yECP (t)〉 = y0(t), que e o valor esperado da ECP.

Nas simulacoes foi considerado y0(0) = 1 e o erro foi calculado em t = 1.



Muito bla bla bla . . .Mas como isso funciona na pratica

Na pratica tem-se a equacao diferencial

dy(t)

dt= −ky(t)

A natureza estocastica da variavel k ≡ k(ω) e conhecida. Por exemplo, ke uma variavel aleatoria Gaussiana com media µk e variancia σ2

k.

Escolhe-se o conjunto de polinomios Φi∞i=0 da ECP.

O resto e conta!

E importante

Escolher o conjunto de polinomios que esteja relacionado com as variaveisaleatorias do problema.

Escolher adequadamente os polinomios:

facilita o calculo dos produtos internos e dos coeficientes;garante convergencia exponencial da solucao (veja nosproximos slides);



Como encontrar os coeficientes da expansao

Os coeficientes da expansao veem do teorema de series de Fourier, so quena serie truncada.

k(ω) =

P∑i=0

〈k(ω), Φi(ξ(ω))〉〈Φi, Φi〉

Φi(ξ(ω))

Encontrar 〈k(ω), Φi(ξ(ω))〉 depende da variavel aleatoria representadapor k. Na pratica, esse produto interno e uma integral.

〈k(ω), Φi(ξ(ω))〉 =

∫D

kΦi(ξ)dP (ξ)

Formas de resolver essa integral:

na racametodos numericos (quadraturas)Monte–Carlo



Exemplo de como encontrar os coeficientes da expansao

Considere k uma variavel aleatoria exponencial, logo, f(k) = e−k parak > 0.

A inversa da sua funcao distribuicao F (k) e

k = h(u) ≡ F−1(u) = − ln(1− u)

onde u ∼ U(0, 1)

Usando ξ ∼ Gaussiana(0, 1) tem-se que a sua inversa e:

ξ = l(u) ≡ G−1(u)

Substituindo na integral do produto interno

〈k(ω), Φi(ξ(ω))〉 =

∫ 1

0

h(u)Φi(l(u))du

E agora e maos a obra!



Distribuicao Gaussiana e polinomios de Hermite

k ∼ Gaussiana(k; 0, 1)

Φi sao polinomios de Hermite.

Figura: Convergencia do erro no valor esperado.



Distribuicao Gama e polinomios de Laguerre

k ∼ Gama(k; α)

Φi sao polinomios de Laguerre.

Figura: Convergencia do erro no valor esperado. α = 0 distribuicao exponencial (quadrados), α = 1(triangulos).



Distribuicao Poisson e polinomios de Charlier

k ∼ Poisson(k; λ)

Φi sao polinomios de Charlier.

Figura: Convergencia do erro no valor esperado. λ = 1 (quadrados), λ = 2 (triangulos).



Distribuicao Gama e polinomios de Hermite

k ∼ Exponencial(k; 1)

Φi sao polinomios de Hermite e de Laguerre.

Figura: Convergencia do erro no valor esperado nao e exponencial. Polinomios de Hermite (quadrados),polinomios de Laguerre α = 0 (triangulos).



Comparacao com metodo de Monte–Carlo

Com a ECP observa-se que o erro e da ordem de 10−3 com valoresde P = 2.

Para P = 4 obtem-se erros da ordem de 10−4 (polinomios deHermite) a 10−9 (polinomios de Jacobi).

Uma simulacao de Monte–Carlo na mesma equacao diferencial queconsidera k ∼ Gama(k; 0) (distribuicao exponencial) apresenta osseguintes resultados:

N 102 103 104 105

ε 4.0× 10−2 1.1× 10−2 5.1× 10−3 6.5× 10−4

Tabela: Convergencia do erro no valor esperado para a simulacao de Monte–Carlo.


Oscilador aleatorio de segunda ordem

Oscilador com excitacoes aleatorias

Considere o sistema de equacoes diferenciaisdx(t)

dt = y(t)dy(t)

dt + cy(t) + kx(t) = f(t)

Assume-se que

c ≡ c(ω) = c + σcξ1(ω)k ≡ k(ω) = k + σkξ2(ω)f(t) ≡ f(t, ω) = (f + σfξ3(ω)) cos(wt)As variaveis aleatorias sao Gaussianas com media 0 e variancia1 e sao independentes

Tem-se portanto que

x(t) ≡ x(t, ω)y(t) ≡ y(t, ω)




A ECP aplica-se ao vetor aleatorio ξ = (ξ1, ξ2, ξ3).

A expressao geral para os polinomios de Hermite e dada por

Hn(ξi1 , . . . , ξin) = e12 ξT ξ(−1)n ∂n

∂ξi1 . . . ∂ξin

e12 ξT ξ

essa representacao tambem e conhecida como formula de Rodriguez

Aplicando a ECP a x(t, ω)

x(t, ω) = x0(t)H0 +3∑

i=1

x1i(t)H1(ξi)

+3∑

i=1

i∑j=1

x2ijH2(ξi, ξj)




Como todos os termos sao ortogonais entre si adota-se uma notacaoreduzida

x(t, ω) =∑

i xi(t)Φi(ξ)= x0(t)H0 + x1(t)H1(ξ1) + x2(t)H1(ξ2) + x3(t)H1(ξ3)+

x4(t)H2(ξ1, ξ1) + x5(t)H2(ξ2, ξ2) + x6(t)H2(ξ3, ξ3)+x7(t)H2(ξ1, ξ2) + x8(t)H2(ξ2, ξ3) + . . .

esta serie e truncada em P termos.

As variavais aleatorias ξ1, ξ2, ξ3, tambem sao representadas nessabase.

c =∑

i ciΦi(ξ) = c + σcH1(ξ1)k =

∑i kiΦi(ξ) = k + σkH1(ξ2)

c = cos(wt)(∑

i fiΦi(ξ)) = cos(wt)(f + σfH1(ξ3))




Aplicando a ECP a equacao diferencial∑

i

dxi(t)

dtΦi =

∑i

yi(t)Φi

∑i

dyi(t)

dtΦi +

∑i

∑j

ciyj(t)ΦiΦj +∑

i

∑j

kixj(t)ΦiΦj =∑

i

fi(t)Φi

Aplicando o metodo de Galerkin tem-se o sistema de equacoesdiferenciais deterministicas

dxi(t)

dt= yi(t)

dyi(t)

dt+

1

〈Φi, Φi〉∑

i

∑j

(ciyj(t) + kixj(t))〈ΦiΦj , Φk〉 = fi(t)

Agora e calcular as integrais de produto interno e rodar a maquina defazer salcicha.



Comparando o erro da ECP com a solucao esperada

Comparando o erro entre o valor esperado da ECP (x0(t)) com o valoresperado da solucao do sistema.

Figura: Convergencia do erro no valor esperado nao e exponencial.



O numero de slides e finito!Obrigado pela paciencia.


[Dongbiu Xiu et al]The Wiener-Askey Polynomial Chaos For StochasticDifferential EquationsSIAM Journal on Scientific Computing, vol. 24, no. 2, 2002

[Dongbiu Xiu et al]Stochastic Modeling of Flow-Structure Interactions usingGeneralized Polynomial ChaosDivision of Applied Mathematics, Brown University,September, 2001

[Andrew J. Newman]Model reduction via the Karhunen-Loeve Expansion. Part I: AnExpositonInstitute for Systems Research and Eletrical EngineeringDepartment, University of Maryland, April, 1996

[Carlos Kubrusly]Elements of Operator Theory


Birkhauser, 2001

[Roger G. Ghanem et al]Stochastic finite elementsDover, 1991

Education

Expansão em caos polinomial