André Fernandes Vasconcelos 21 de Julho de 2007calves/AN/lmac/Projecto/... · bem como aperfeiçoar os métodos numéricos que nos permitem prever a na- tureza, localização, e/ou

Detecção de cavidades ou inclusões através damedição de condutividade na superfície

exterior -O Problema Directo

André Fernandes Vasconcelos

21 de Julho de 2007

1

Conteúdo1 Introdução 3

2 Formulação e Teoria 42.1 Formulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.1.1 Problema Directo . . . . . . . . . . . . . . . . . . . . . 42.1.2 Problema Inverso . . . . . . . . . . . . . . . . . . . . . 4

2.2 Existência e Unicidade . . . . . . . . . . . . . . . . . . . . . . 52.3 Resultados de Análise Funcional . . . . . . . . . . . . . . . . . 62.4 Densidade . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72.5 O Método das Soluções Fundamentais . . . . . . . . . . . . . . 11

2.5.1 Solução aproximada . . . . . . . . . . . . . . . . . . . 112.5.2 Escolha dos pontos fonte . . . . . . . . . . . . . . . . . 122.5.3 Métodos dos Mínimos Quadrados . . . . . . . . . . . . 132.5.4 Regularização de Tikhonov . . . . . . . . . . . . . . . . 14

3 Simulações 153.1 Concretização . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

3.1.1 O Problema . . . . . . . . . . . . . . . . . . . . . . . . 153.1.2 O Método das Soluções Fundamentais . . . . . . . . . 153.1.3 Ruído e Erro . . . . . . . . . . . . . . . . . . . . . . . 163.1.4 Introdução da fronteira Γ∗2 . . . . . . . . . . . . . . . . 183.1.5 Matlab . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Resultados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2.1 Previsões . . . . . . . . . . . . . . . . . . . . . . . . . 203.2.2 Domínio . . . . . . . . . . . . . . . . . . . . . . . . . . 203.2.3 Proximidade das fronteiras auxiliares . . . . . . . . . . 223.2.4 Introdução da fronteira auxiliar Γ∗2 . . . . . . . . . . . 243.2.5 Ruído . . . . . . . . . . . . . . . . . . . . . . . . . . . 253.2.6 Oscilação . . . . . . . . . . . . . . . . . . . . . . . . . 27

4 Conclusões 28

5 Anexo 305.1 ProblemaDirecto.m . . . . . . . . . . . . . . . . . . . . . . . . 30

2

1 IntroduçãoA capacidade de detectar danos internos ou anomalias numa estrutura

sem a destruir ou afectar de qualquer forma é de uma extrema importânciapara campos como os da Medicina e da Engenharia Civil. Nesse sentido, éimportante não só aperfeiçoar o equipamento que se utiliza para ler dados,bem como aperfeiçoar os métodos numéricos que nos permitem prever a na-tureza, localização, e/ou forma das anomalias a partir dos dados fornecidos.

Esse é o tema que abordamos no ramo de Problemas Inversos para asEquações Diferenciais Parciais; no nosso caso, o Problema Inverso (a razãodeste nome será evidente mais tarde) será a detecção de uma cavidade (porexemplo, numa viga), medindo a condutividade (de calor) na superfície ex-terna. No entanto, para abordarmos o Problema Inverso deveríamos primeiroter um entendimento completo do Problema Directo, que consiste em, dadaa cavidade, descobrirmos a condutividade gerada.

Neste Projecto estudamos o Problema Directo. No capítulo ‘Formulação eTeoria’ definiremos o Problema Directo, o Problema Inverso, e definiremos osmétodos numéricos utilizados para resolver o Problema Directo, estando entreestes o Método das Soluções Fundamentais, e a Regularização de Tikhonov;ainda provamos um resultado de densidade que prova a validade do nossométodo. Finalmente, no capítulo ’Simulações’ pomos em prática tudo o quefoi introduzido e são apresentados os resultados, que são depois discutidosnas Conclusões.

3

2 Formulação e Teoria

2.1 Formulação

2.1.1 Problema Directo

O nosso Problema Directo consiste no seguinte: Temos um domínio Ω (biou tridimensional), que representa o nosso meio homogéneo, e um domínio ω,com ω ∈ Ω, que representa a cavidade no nosso domínio. O nosso problemaé determinar u tal que

∆u = F, Ω\ωu = g, γ = ∂ω

u = G, Γ = ∂Ω

onde F ∈ L2(Ω), g ∈ C(γ), G ∈ C(Γ), são funções dadas. De agora emdiante, considera-se F = 0.

2.1.2 Problema Inverso

O Problema Inverso é formulado da seguinte forma: seja u tal que

∆u = F, Ω\ωu = c, c ∈ R

e tal que conhecemos ∂nu(x) = g1(x) sobre parte de Γ, e temos o domínioΩ (bi ou tridimensional) conhecido. Temos como incógnita ω (a sua local-ização e forma), e o nosso problema é determiná-lo. No nosso Projecto nãoiremos examinar este problema, apresentando-o apenas porque o Projectoestá inserido no contexto deste Problema.

4

2.2 Existência e Unicidade

É necessário saber se temos condições para que a solução exista e sejaúnica, no nosso Problema Directo. Vemos que, segundo a teoria clássica dasEDP’s, é esse o caso:

Teorema 1 Consideremos o Problema Directo definido acima. Se Γ e γforem C1, então a solução clássica para o Problema Directo existe, é única,e pertence a C(Ω) ∩ C2(Ω).

Temos ainda que:

Teorema 2 Suponhamos que estamos nas condições do Teorema 1, e aindaque g ∈ C∞(γ), G ∈ C∞(Γ). Então a solução clássica para o ProblemaDirecto existe, é única, e é de classe C∞(Ω).

No entanto a teoria clássica de EDP’s pode ser muito restritiva para onosso contexto. Incluímos no grupo das nossas soluções as soluções fracas(ver [2]), e obtemos o seguinte resultado:

Teorema 3 Suponhamos que g ∈ H1/2(γ) e G ∈ H1/2(Γ). Então a soluçãofraca ∈ H1(Ω\ω) existe, e é única.

Introduzimos ainda a Solução Fundamental da equação de Laplace (∆u =0). Esta é dada por Φ(x) = 1

2πlog |x|, se estivermos em 2 dimensões, e

Φ(x) = 14π|x| , se estivermos em 3 dimensões (as Soluções Fundamentais são

tratadas em detalhe em [2]).

5

2.3 Resultados de Análise Funcional

Na secção seguinte, na qual introduziremos o Método das Soluções Fun-damentais, necessitaremos de algumas noções de Análise Funcional, bemcomo do Teorema da Continuação Analítica. Introduzimo-los a todos aqui,podendo esta subsecção ser lida apenas como referência para a próxima sub-secção.

Definição 1 Sejam X, Y espaços de Hilbert e T : X → Y um operadorlinear limitado. Então o operador adjunto T ∗ de T é o operador T ∗ : Y → Xtal que, para quaisquer x ∈ X e y ∈ Y ,

< Tx, y >X=< x, T ∗y >Y

onde < ., . >X e < ., . >Y são os produtos internos definidos em X e de Y ,respectivamente. Regra geral omite-se o X e o Y no produto interno, sendoeles identificados pelo contexto.

Definição 2 Sejam X um espaço de Hilbert e A um seu subconjunto nãovazio. Designamos por ortogonal de A o conjunto A⊥ formado pelos elemen-tos de X que são ortogonais a todos os elementos de A, i.e.

A⊥ = x ∈ X : x ⊥ A

Proposição 4 ker T ∗ = (r(T ))⊥

Demonstração:Dados x ∈ X e y ∈ ker T ∗ temos

0 =< x, T ∗y >=< Tx, y >

pelo que ker T ∗ ⊂ (r(T ))⊥ (porque qualquer elemento de ker T ∗ é ortogonala Tx, x ∈ X). Por outro lado, se y ∈ (r(T ))⊥, e x = T ∗y, tem-se

0 =< TT ∗y, y >=< T ∗y, T ∗y >= ‖T ∗y‖,

o que significa que y ∈ ker T ∗, o que conclui a nossa demonstração.

2

Teorema 5 Dado um subespaço vectorial A 6= ∅ de um espaço de Hilbert X,A é denso em X sse A⊥ = 0.

Acabamos a secção com o seguinte Teorema, útil para provar um resultadode densidade:

Teorema 6 (Teorema de Continuação Analítica) Seja u uma funçãoanalítica em W ⊃ ω (ω aberto). Os valores de u em ω determinam u em Wde forma única.

6

2.4 Densidade

Nesta secção provamos um resultado acerca da densidade das SoluçõesFundamentais (definidas no final da subsecção 2.2). Esta subsecção podeser lida após a introdução do Método das Soluções Fundamentais, na secçãoseguinte. O nosso resultado usa o seguinte lema:

Lema 7 Sejam Ω, Γ, γ definidos como no Problema Directo, e ainda umaberto U ⊂ Rn\(Ω\ω). Definimos o operador A : L2(∂U) → L2(Γ ∪ γ),tal que

Af(x) =

∫∂U

f(y)Φ(x− y)dy.

Então A∗g(y) =∫

Γ∪γg(x)Φ(x− y)dx.

Demonstração:O nosso objectivo é encontrar A∗ tal que < Af, g >L2(Γ∪γ)=< f, A∗g >L2(∂U),

f ∈ L2(∂U), g ∈ L2(Γ ∪ γ). Temos que

< Af, g >L2(Γ∪γ) =

∫Γ∪γ

g(x)

∫∂U

f(y)Φ(x− y)dydx

=

∫∂U

f(y)

∫Γ∪γ

Φ(x− y)g(x)dxdy

=< f,

∫Γ∪γ

Φ(x− y)g(x)dx >L2(∂U)

2

Note-se que spanΦ(x− y)|Γ∪γ, y ∈ ∂U

= r(A). O resultado de densidade

para 3D é dado pelo seguinte teorema:

Teorema 8 Considere-se o Problema Directo formulado, e consideremosum aberto U tal que U ⊂ Rn\(Ω\ω) (n = 3), tal que U tem pelo menosum ponto em cada componente conexa de Rn\(Ω\ω). Então L2(Γ ∪ γ) =

spanΦ(x− y)|Γ∪γ, y ∈ ∂U

.

Demonstração:Queremos provar que r(A) é denso em L2(Γ ∪ γ). Segundo o Teorema 5,

para isso temos que provar que r(A)⊥ = 0, o que, segundo o Proposição 4é equivalente a provar ker A∗ = 0. Recordemos que temos, pelo nosso Lema,que A∗g(y) =

∫Γ∪γ

g(x)Φ(x− y)dx.Então temos a provar que, se ug := A∗g ≡ 0 em Γ∗ ∪ γ∗, temos que

ug ≡ 0 em todo o R2, o que implica g ≡ 0. A conclusão desta demonstraçãoé semelhante à do teorema seguinte.

7

2

Se este teorema se verificasse em 2D, saberíamos que basta escolhermos pon-tos yi em curvas envolvendo o nosso domínio Ω\ω.

No entanto, este teorema só se aplica no caso 3D, porque a solução fun-damental nessa dimensão decresce para zero no infinito. No caso 2D, a nossasolução fundamental tem crescimento logarítmico, e há pelo menos um passoda demonstração do teorema anterior que falha; necessitamos de enfraquecero resultado.

Teorema 9 Considere-se o Problema Directo formulado, e consideremosum aberto U tal que U ∈ Rn\(Ω\ω) (n = 2), tal que U tem pelo menosum ponto em cada componente conexa de Rn\(Ω\ω). Então [L2(Γ ∪ γ)] =

spanΦ(x− y)|Γ∪γ, y ∈ ∂U

, onde [L2(Γ ∪ γ)] é o espaço L2 onde introduz-

imos a relação de equivalência f ≡ g sse f = g + C, C ∈ R.

Temos que alterar o produto interno em função do novo espaço introduzido,e em consequência teremos também que alterar A∗.

Proposição 10 O funcional bilinear definido por

< [f ], [g] >=

∫Γ∪γ

f ∗(x)g∗(x)dx,

com f ∗ o representante de [f ] tal que∫

Γ∪γf ∗(x)dx = 0, define um produto

interno em [L2](Γ ∪ γ).

Note-se que o representante de f ∗ de [f ] está bem definido, bem como oproduto interno. Sendo assim, o funcional A é mantido inalterado, e A∗

toma a seguinte forma:

A∗[g](y) =

∫Γ∪γ

g∗(x)Φ(x− y)dx

. Demonstração (Teorema 9):Semelhantemente ao Teorema 8, temos a provar que, se u[g] := A∗[g] ≡ 0

em Γ∗ ∪ γ∗, temos que u[g] ≡ 0 em todo o R2, o que implica [g] ≡ 0. Sejaentão [g] tal que A∗[g](x) =

∫Γ∪γ

Φ(x− y)g∗(x)dx ≡ 0. A nossa resoluçãodepende do aberto U escolhido: supomos sem perda de generalidade que U écomposto pelos abertos A e E (representados abaixo), sendo a sua fronteiracomposta por Γ∗ e γ∗. Note-se que Ω = C ∪D ∪ E, ω = D ∪ E.

8

Figura 1: O domínio considerado para a nossa demonstração

i. Comecemos por provar que u[g] é nulo em E. Temos que u[g] obedece àseguinte equação:

∆u[g] = 0, E

u[g] = 0, ∂E

Então pelo Teorema 1 existe apenas uma solução, que é dada por u[g] ≡0 em E.

ii. Provamos que u[g] é nulo em D. u[g] é analítica em R2\(Γ∪ γ), e esta énula em E; logo temos pelo Teorema da Continuação Analítica que sepode estender u[g] ≡ 0 para D.

iii. Provamos que u[g] é nulo em A. Temos que u[g] obedece à seguinteequação:

∆u[g] = 0, A

u[g] = 0, Γ∗

u[g](x) = α log(|x|) + O(1), α =

∫γ∪Γ

g∗(y)dy, |x| → ∞

Para cada α fixo (é este o caso), temos que este problema tem soluçãoúnica. No nosso caso, e devido a todo o trabalho que tivemos a definirA∗, temos que α = 0. Então temos u[g] ≡ 0 como solução única destaequação.

9

iv. Provamos que u[g] é nulo em B. Mais uma vez utilizamos o Teoremada Continuação Analítica, para, a partir de A, “estendermos” a funçãonula u[g] ≡ 0 para B.

v. Resta-nos provar que u[g] é nulo em C. Repare-se que do lado de “fora”de C, temos u[g] ≡ 0 em γ e Γ; chamamos à solução exterior u−[g] (àinterior u+

[g]). u[g] é o que se chama um Potencial de Camada Simples,e existem resultados que garantem que este é contínuo através de γ eΓ (ver [fonte]). Logo u+

[g] ≡ 0, e temos que

∆u[g] = 0, C

u[g] = 0, Γ ∪ γ

o que evidentemente só tem solução nula u[g] ≡ 0.

Logo temos u[g] ≡ 0, em R2\(Γ ∪ γ), o que significa que g ≡ 0 e portanto[g] ≡ 0. Logo ker A∗ = 0.

2

10

2.5 O Método das Soluções Fundamentais

2.5.1 Solução aproximada

Recordemos que o nosso problema é dado por

∆u = 0, Ω\ωu = g, γ = ∂ω

u = G, Γ = ∂Ω,

em que vamos definir ω e Ω como sendo elipses (tal como na Figura 1).Nós sabemos que, se y /∈ Ω\ω, então uy(x) := Φ(x− y) é uma solução da

equação de Laplace. Então será natural pensar que poderemos aproximar anossa solução por u(x) =

∑ni=1 αiuyi

(x) =∑n

i=1 αiΦ(x− yi), com yi /∈ Ω\ω(aos quais chamaremos pontos fonte) - restaria escolher as constantes αi demodo a escolhermos a melhor aproximação da solução exacta na fronteira.Este raciocínio está de acordo com os resultados de densidade da secçãopassada, pois span

Φ(x− y)|Γ∪γ, y ∈ ∂U

= r(A) (com U definido como no

Teorema 8).Para ser preciso, pretendemos minimizar o erro da nossa aproximação na

norma L2(Γ ∪ γ): esta surge como a maneira natural de medir o erro, nocontexto dos resultados de densidade dados anteriormente, e esta reflecte aqualidade da aproximação da solução exacta no domínio inteiro; pelo Princí-pio do Máximo, para a equação de Laplace o máximo do erro tem que estarna fronteira.

Para minimizar o erro da solução aproximada na fronteira, na normaL2(Γ ∪ γ), recorremos ao Método dos Mínimos Quadrados.

11

2.5.2 Escolha dos pontos fonte

Uma questão crucial para a formulação do nosso Método dos MínimosQuadrados é a seleção dos pontos yi, e das fronteiras Γ∗ e γ∗ (ver a sub-secção 2.4) sobre as quais se situam esses pontos. Necessitamos de escolherU para sabermos a forma das nossas fronteiras Γ∗ e γ∗, tal que estejamos nascondições do Teorema 9.

A nossa primeira escolha é escolher Γ∗ e γ∗ como na Figura 1; embora nãotenhamos um resultado de densidade completa, consideramos que o facto determos densidade a menos de constante é tão bom como ter densidade.

Se escolhermos Γ∗ = Γ∗1 ∪ Γ∗2, em que Γ∗1 é um anel que envolve Γ, eΓ∗2 é um anel que envolve Γ∗1, conseguimos um resultado de densidade emL2(Γ∪ γ), o que poderá significar melhores resultados para o mesmo númerode pontos fonte utilizados.

É de notar que não temos quaisquer resultados sobre a convergência donosso método em função do número de pontos fonte utilizados. Veremosmais tarde que, mesmo havendo qualquer convergência, esta seria anuladapor erros de mau condicionamento.

12

2.5.3 Métodos dos Mínimos Quadrados

Escolhemos então nx +nz valores conhecidos (digamos, g(x1), . . . , g(xnx),G(z1), . . . , G(znz)), tão equidistribuídos quanto os domínios nos permitam.Utilizamos o MMQ para determinar as constantes αi, i = 1, . . . , n, tal que oerro da nossa solução aproximada u na norma discreta l2 seja mínimo; paraisso temos que deduzir o sistema linear correspondente. Suponhamos queera possível interpolarmos os nx + nz pontos com valores conhecidos. Entãoteríamos o seguinte sistema linear sobredeterminado:

Φ(x1 − y1) . . . Φ(x1 − yn)...

...Φ(xnx − y1) . . . Φ(xnx − yn)Φ(z1 − y1) . . . Φ(z1 − yn)

......

Φ(znz − y1) . . . Φ(znz − yn)

α1

...αn

=

g(x1)...

g(xnx)G(z1)

...G(znz)

No entanto o sistema é, geralmente, impossível (especialmente se tiver-

mos introduzido ruído nos dados). Chamando à matriz M , o nosso MMQ éaplicado ao multiplicarmos o sistema por MT à esquerda. Desta forma obte-mos um sistema possível (as funções Φ(x− yi), i = 1, . . . , n são linearmenteindependentes: ver a página 130 de [1]) com matriz quadrada MT M , e quecorresponde a uma minimização das distâncias na norma l2. Temos entãoque resolver o sistema

MT M

α1...

αn

= MT

g(x1)...

g(xnx)G(z1)

...G(znz)

No entanto o sistema resultante poderá ser mal condicionado, gerando errosnuméricos relevantes. Para isso utilizamos a Regularização de Tikhonov.

13

2.5.4 Regularização de Tikhonov

A Regularização de Tikhonov está inserida no contexto dos EsquemasRegularizadores, que procuram aproximar a inversa de um operador injectivoe compacto, mesmo que geralmente este seja não sobrejectivo e portanto nãotenha uma inversa definida em todo o espaço de chegada.

No nosso caso o nosso operador é o operador linear M : Rn → Rnx+nz .A Regularização de Tikhonov aplicada ao nosso sistema linear é definida daseguinte forma:

(αI + MT M)

α1...

αn

= MT

g(x1)...

g(xnx)G(z1)

...G(znz)

, α > 0

A Regularização de Tikhonov reduz-se ao MMQ quando escolhemos α =0, i.e., reduz-se à solução exacta (caso esta exista!) que pretendíamos.

Encontrar um α óptimo é uma tarefa não trivial. Por um lado, um αmuito pequeno deixa o sistema mal condicionado, e, caso tenhamos intro-duzido ruído nos dados da fronteira, este estraga a nossa solução; mas afidelidade à solução exacta é mantida. Por outro lado, um α grande tornao sistema bem condicionado e diminui os efeitos do ruído, mas afastamo-nosda solução que procuramos.

Existe um resultado que nos garante que existe um resultado óptimo:

Teorema 11 (Princípio da Discrepância aplicado a Tikhonov) Seja A :X → Y um operador linear, injectivo, compacto, e A(X) = Y . Seja g ∈A(X) e gδ ∈ Y tais que ‖g− gδ‖ ≤ δ < ‖gδ‖. Então existe um único α(δ) talque ‖ARα(δ)gδ − gδ‖ ≤ δ e temos Rα(δ)gδ → A−1gδ, quando δ → 0.

No entanto, não existe nenhuma maneira determinística de escolher o α óp-timo; no nosso caso fazemos experimentação não exaustiva até obter umα tão pequeno quanto possível que nos dê um sistema bem condicionado;entendemos por bem condicionado um sistema linear que tenha número decondição de ordem superior a 108.

14

3 Simulações

3.1 Concretização

Uma vez definido o método, é preciso definir mais concretamente o prob-lema; de seguida iremos implementar o algoritmo.

3.1.1 O Problema

Primeiro que tudo é necessário definir o problema:

• Os abertos Ω e ω;

• As funções g e G.

Os abertos Ω e ω serão dados por elipses, de eixos variáveis, respectivamente,(Rx, Ry) e (rx, ry). Iremos testar se o algoritmo sofre diferenças significativascom os diferentes domínios, e depois fixaremos um domínio para estudaroutras variáveis.

Por outro lado g e G serão gerados de duas maneiras: utilizando soluçõesexactas harmónicas já conhecidas, ou utilizando qualquer função (o que sig-nifica que não poderemos comparar visualmente as soluções exactas com asaproximadas, em todo o domínio). No primeiro grupo das funções temos

• h1(x, y) = x2 − y2

• h2(x, y) = log(x2 + y2)

No segundo grupo temos

• i1(r, θ) = 1 + r cos(θ)

• i2(r, θ) = 1 + r cos(2θ)

• i3(r, θ) = 1 + r cos(3θ)

3.1.2 O Método das Soluções Fundamentais

Precisamos de definir:

• As fronteiras auxiliares γ∗ (elipse de eixos rx e ry) e Γ∗ (elipse de eixosRx e Ry);

• O número de pontos nx e nz utilizados em γ e Γ, respectivamente;

• O número de pontos n1 e n2 utilizados em γ∗ e Γ∗, respectivamente;

15

As fronteiras auxiliares γ∗ e Γ∗ irão ser, respectivamente, proporcionais a γ(por um factor q1) e Γ (por um factor q2). Escolheremos n1 = n2 = n, e nx =nz = 2n, variando o n, e a a distribuição de pontos ao longo de cada fronteiraserá equidistante. Por exemplo, temos xi = (rx cos(2(i−1)π/nx), ry sin(2(i−1)π/nx)), e zi = (Rx cos(2(i−1)π/nz), Ry sin(2(i−1)π/nz)), com i não nec-essariamente inteiro. O caso será ligeiramente diferente quando adicionarmosmais fronteiras auxiliares.

3.1.3 Ruído e Erro

Iremos também introduzir ruído nos nossos valores lidos. O ruído intro-duzido será proporcional, por um factor ε, a ‖G‖2

L2(Γ), para Γ, e ‖g‖2L2(γ),

para γ.O erro relativo irá ser dado por

‖u− u‖L2(Γ∪γ)

‖u‖L2(Γ∪γ)

=

√√√√∫Γ∪γ

(u(x)− u(x))2dx∫Γ∪γ

u(x)2dx

que por sua vez também terá que ser aproximado numericamente. Poderíamosser tentados a aproximar

∫Γ∪γ

(u(x)− u(x))2dx por

n∑i=1

ci

2((u(xi+1)− u(xi+1))

2 + (u(xi)− u(xi))2)

+n∑

i=1

Ci

2((u(zi+1)− u(zi+1))

2 + (u(zi)− u(zi))2)

onde xi e zi são os pontos com valores conhecidos, e ci = |xi+1 − xi| eCi = |zi+1 − zi|; corresponde à utilização da regra dos trapézios numa curvaaproximada por segmentos de recta. No entanto esta aproximação não é cor-recta, pois estamos a aproximar o integral do erro usando exactamente ospontos onde o tentamos minimizar (pois estamos a utilizar Mínimos Quadra-dos, ou a Regularização de Tikhonov). Uma melhor aproximação será dadapor

n∑i=1

ci+1/2

2((u(xi+1/2)− u(xi+1/2))

2 + (u(xi−1/2)− u(xi−1/2))2)

+n∑

i=1

Ci+1/2

2((u(zi+1/2)− u(zi+1/2))

2 + (u(zi−1/2)− u(zi−1/2))2)

16

onde ci+1/2 = |xi+1/2 − xi−1/2| e Ci+1/2 = |zi+1/2 − zi−1/2|. Esta continua aser uma aproximação de primeira ordem do integral.

A aproximação para∫

Γ∪γu(x)2dx é semelhante à anterior.

17

3.1.4 Introdução da fronteira Γ∗2

No capítulo sobre Densidade, ficou sugerido que melhores resultados pode-riam ser obtidos se introduzíssemos uma terceira fronteira auxiliar Γ∗2, envol-vendo Γ∗2. Criamo-la proporcional a Γ∗, por um factor de 3/2. Os pontosfonte colocados em Γ∗2 estarão colocados com um ângulo adicional de π/n2(em que n2 é o número de pontos em Γ∗, e em Γ∗2). A diferença entre as duasabordagens é representada graficamente no exemplo da figura abaixo.

Figura 2: Temos Γ a azul, Γ∗ a vermelho (com 20 pontos fonte, ’+’), ea segunda abordagem com Γ∗ e Γ∗2 a verde (com 20 pontos fonte ’o’). Afronteira γ não está representada.

18

3.1.5 Matlab

A implementação em Matlab de tudo o que estivemos a preparar atéesta secção é realizada na função ProblemaDirecto, com auxílio das funçõesgeragG, e funcaoteste. A função recebe os argumentos (rx, ry, q1, n1, g, Rx,Ry, q2, n2, G, ruido, alpha, anel):

• (rx,ry) são os comprimentos do eixo de γ, e (Rx,Ry) são os comprimen-tos do eixo de Γ;

• ’q1’ e ’q2’ são as proporções respectivas de γ∗ em relação a γ, Γ∗ emrelação a Γ.

• ’n1’ e ’n2’ são o número de pontos em γ∗ e Γ∗, respectivamente; são ospontos yi;

• ’g’ e ’G’ são as listas de (’nx’ e ’ny’) valores em γ e Γ, respectivamente.São construídas pela função geragG, que deverá ser corrida antes dafunção ProblemaDirecto.

• ’ruido’ é a proporção de ruído introduzida nos dados iniciais;

• ’alpha’ é o parâmetro de regularização de Tikhonov;

• ’anel’ define se utilizamos uma segunda fronteira auxiliar (Γ∗2), caso noqual valerá 1, ou se não a utilizamos, caso no qual valerá 0.

A partir desses dados, a função ProblemaDirecto devolve os valores da soluçãoaproximada nos pontos xi+1/2, zi+1/2, i = 1, . . . , n, da solução exacta nos mes-mos pontos, a lista de valores de αi, e o erro relativo da solução aproximada.Para ver mais detalhes da função, dever-se-á consultar o código comentadono Anexo.

19

3.2 Resultados

3.2.1 Previsões

Esperamos obter os seguintes resultados:

• O erro deverá diminuir com o número de pontos, a menos de maucondicionamento. A velocidade de decrescimento é desconhecida;

• Variações do domínio não afectarão demasiado a eficiência do algo-ritmo;

• A proximidade das fronteiras auxiliares deverá afectar o resultado, emtermos de erro e mau condicionamento;

• O algoritmo deverá ser bastante sensível ao ruído;

• Funções com grandes oscilações na fronteira deverão piorar os resulta-dos;

• A introdução da fronteira Γ∗2 deverá dar melhores resultados.

3.2.2 Domínio

Testamos os domínios com as seguintes características:

i. (rx, ry) = (2, 1), (Rx, Ry) = (4, 4);

ii. (rx, ry) = (2, 2), (Rx, Ry) = (4, 4);

iii. (rx, ry) = (2, 1), (Rx, Ry) = (6, 6);

iv. (rx, ry) = (3, 2), (Rx, Ry) = (4, 4);

O primeiro exemplo será o nosso exemplo padrão; o segundo exemplo testase concentricidade dos domínios afectará os resultados, o terceiro e o quartotestam se a distância entre domínios afecta os resultados.

Recorde-se que nx = nz = 2n, n1 = n2 = n. Utilizamos q1 = 1/8, q2= 4, ruido = 0, e alpha variável. Os resultados são apresentados na seguintetabela de erros relativos. Marcamos soluções resultantes de sistemas malcondicionados, e que utilizaram portanto a Regularização de Tikhonov, comum asterisco (*).

20

h1(x, y) = x2 − y2

Domínio \ n 5 10 20 401 1.6324e-001 6.0501e-005 4.5829e-008* 3.3080e-008*2 1.6199e-001 5.9947e-005 2.2874e-008* 4.9989e-008*3 1.6402e-001 6.0826e-005 9.7172e-008* 2.5831e-008*4 1.5903e-001 5.7598e-005 2.8989e-007 5.0824e-008*

h2(x, y) = log(x2 + y2)Domínio \ n 5 10 20 40

1 1.6458e-004 3.0388e-007 2.2657e-009* 1.3877e-009*2 9.1427e-011 5.4182e-012 2.7761e-010* 7.3562e-010*3 4.5561e-003 1.4519e-004 1.1784e-007* 1.1522e-009*4 3.7945e-003 9.8468e-005 3.4963e-008 1.0542e-009*

i1(r, θ) = 1 + r cos(θ)Domínio \ n 5 10 20 40

1 3.3442e-003 1.4365e-006 4.3518e-009* 3.6235e-009*2 3.2709e-003 1.4042e-006 5.2553e-009* 4.4106e-009*3 3.5793e-003 1.5382e-006 1.2698e-008* 1.2085e-008*4 3.0723e-003 1.3037e-006 2.2990e-008 4.5366e-009*

Examinando as tabelas, chegamos à conclusão de que, com uma excepção,os resultados são qualitativamente invariantes pelo domínio, embora umamenor distância entre Γ e Γ∗ torne as matrizes mais bem condicionadas (nocaso 4). A única excepção ocorre no caso em que os domínios são concêntricos(caso 2), o que proporciona excelentes resultados para uma função radial(h2). Nota-se também que a convergência em função do número de pontos éestragada pelo mau condicionamento.

21

3.2.3 Proximidade das fronteiras auxiliares

Fixamos (rx, ry) = (2, 1), (Rx, Ry) = (4, 4) como exemplo canónico, eintroduzimos ruído nulo. Vimos que as distâncias entre fronteira interior eexterior podem provocar mau condicionamento; o mesmo se poderá passarcom as fronteiras auxiliares. Testamos vários pares de q1 e q2.

h1(x, y) = x2 − y2

(q1, q2) \ n 5 10 20 40(1/32, 16) 4.1414e-002 3.6544e-007* 1.2644e-006* 4.3600e-007*(1/16, 8) 8.2601e-002 1.4136e-006 3.0348e-007* 8.7480e-008*(1/8, 4) 1.6324e-001 6.0501e-005 4.5829e-008* 3.3080e-008*(1/4, 2) 3.0809e-001 3.7204e-003 1.5988e-006 1.0710e-008*

(1/2, 3/2) 3.8913e-001 1.8950e-002 1.4095e-004 1.9722e-008(2/3, 4/3) 4.3489e-001 3.5175e-002 8.2026e-004 1.1866e-006(3/4, 5/4) 4.7087e-001 4.9089e-002 2.0674e-003 1.0684e-005(7/8, 9/8) 5.6601e-001 1.0241e-001 8.7110e-003 3.2809e-004

h2(x, y) = log(x2 + y2)(q1, q2) \ n 5 10 20 40(1/32, 16) 2.2436e-006 1.0005e-008* 2.8822e-008* 1.9021e-008*(1/16, 8) 1.8568e-005 4.6173e-009 3.1752e-008* 9.3939e-009*(1/8, 4) 1.6458e-004 3.0388e-007 3.6549e-008* 5.9002e-008*(1/4, 2) 1.5727e-003 2.1213e-005 1.3764e-009 7.1665e-008*


i1(r, θ) = 1 + r cos(θ)(q1, q2) \ n 5 10 20 40(1/32, 16) 5.4383e-005 1.4172e-008* 2.2572e-008* 7.5409e-009*(1/16, 8) 4.3166e-004 7.4017e-009 1.2436e-008* 3.3641e-009*(1/8, 4) 3.3442e-003 1.4365e-006 4.3518e-009* 3.6235e-009*(1/4, 2) 2.3307e-002 3.0829e-004 1.3873e-007 2.4474e-009*


22

Figura 3: Evolução do logaritmo do erro relativo em função de q1, para h1

(verde), h2 (azul), i1 (vermelho), com n = 20. Valores mais baixos significammenor erro relativo.

Os resultados mais consistentes parecem pertencer aos pares (q1, q2) =(1/16,8), (1/8, 4), e (1/4,2), mesmo tendo sistemas mal condicionados. Parafronteiras mais próximas deixamos de ter sistemas mal condicionados, masos erros relativos são piores (os valores das funções de base são muito altosnos pontos zi e xi); para fronteiras mais longínquas o mau condicionamentocomeça a estragar os resultados. Escolheremos como valores padrão para orestante trabalho (q1, q2) = (1/8, 4).

23

3.2.4 Introdução da fronteira auxiliar Γ∗2

Continuamos com (rx, ry) = (2, 1), (Rx, Ry) = (4, 4), e fixamos (q1,q2)= (1/8, 4). Testamos a influência da introdução da fronteira auxiliar Γ∗2nos resultados. Note-se que no caso em que utilizamos a fronteira Γ∗2, temosn1 = n, e n2 = n/2 (arredondado para cima).

h1(x, y) = x2 − y2

Γ∗2 \ n 5 10 20 40Não 1.6324e-001 6.0501e-005 4.5829e-008* 3.3080e-008*Sim 2.4163e-002 9.3116e-005* 1.7550e-005* 2.3432e-008*

h2(x, y) = log(x2 + y2)Γ∗2 \ n 5 10 20 40Não 1.6458e-004 3.0388e-007 2.2657e-009* 1.3877e-009*Sim 1.6578e-004 3.0387e-007* 2.2970e-009* 5.8547e-009*

i1(r, θ) = 1 + r cos(θ)Γ∗2 \ n 5 10 20 40Não 3.3442e-003 1.4365e-006 4.3518e-009* 3.6235e-009*Sim 7.7690e-004 1.8447e-005* 9.0224e-007* 1.7813e-009*

Figura 4: Comparação (do logaritmo) dos resultados em função de n, semΓ∗2 (vermelho), e com Γ∗2 (verde)

A introdução de Γ∗2 chega a dar melhores resultados para alguns casos,mas regra geral o pior condicionamento do sistema faz com que a introduçãoda segunda fronteira auxiliar dê resultados inferiores. Confirma-se que ésuficiente uma fronteira auxiliar.

24

3.2.5 Ruído

No contexto do Problema Inverso, muitas vezes temos de lidar com erroe ruído nos dados introduzidos (pois está associado a leituras de dados emaplicações). Devido a esse facto, consideramos que estudar o efeito de ruídono nosso método para o Problema Directo pode ser útil.

Para este efeito utilizamos (rx, ry) = (2, 1), (Rx, Ry) = (4, 4), (q1, q2)= (1/8, 4).

h1(x, y) = x2 − y2

Ruído \ n 5 10 20 400% 1.6324e-001 6.0501e-005 4.5829e-008* 3.3080e-008*5% 1.7779e-001 8.7092e-002 7.5944e-002* 6.3250e-002*10% 3.0744e-001 1.9938e-001 2.1415e-001* 1.2847e-001*25% 4.9256e-001 5.1212e-001 3.3830e-001* 3.8702e-001*50% 1.1083e+000 1.2439e+000 1.0062e+000* 4.9568e-001*

h2(x, y) = log(x2 + y2)Ruído \ n 5 10 20 40

0% 1.6458e-004 3.0388e-007 2.2657e-009* 1.3877e-009*5% 1.2121e-001 9.0003e-002 1.0437e-001* 8.4144e-002*10% 2.7779e-001 1.9373e-001 1.6319e-001* 8.8492e-002*25% 5.3958e-001 3.3072e-001 3.2361e-001* 2.7017e-001*50% 1.0492e+000 1.1189e+000 9.7539e-001* 8.1795e-001*

i1(r, θ) = 1 + r cos(θ)Ruído \ n 5 10 20 40

0% 3.3442e-003 1.4365e-006 4.3518e-009* 3.6235e-009*5% 1.0738e-001 1.2923e-001 9.7343e-002* 4.9166e-002*10% 2.6581e-001 3.1264e-001 1.5064e-001* 1.3792e-001*25% 6.0225e-001 4.8895e-001 4.7257e-001* 2.6347e-001*50% 1.1796e+000 9.3189e-001 9.3325e-001* 6.8799e-001*

Analisando as tabelas, vê-se que mesmo pequeno ruído reduz drasticamentea qualidade da aproximação. Uma vez introduzido o ruído, o erro aumentaproporcionalmente com o ruído.

25

Figura 5: Evolução do erro relativo em função do ruído, para h1 (verde), h2

(azul), i1 (vermelho), com n = 20.

Figura 6: Soluções aproximada (vermelho) e exacta (azul), para h1, 10% deruído

26

3.2.6 Oscilação

Testamos o nosso método para funções radiais de oscilação crescente.

Ruído \ n 5 10 20 40i1 3.3442e-003 1.4365e-006 4.3518e-009* 3.6235e-009*i2 1.8238e-001 2.8844e-002 1.8742e-002* 1.9044e-002*i3 9.1462e-001 2.6751e-001 5.4497e-002* 4.6487e-002*

Figura 7: Soluções aproximada (ver-melho) e exacta (azul), para i2




As oscilações pioram a qualidade do algoritmo muito rapidamente. Intuiti-vamente, tal deve-se a termos “menos” informação sobre a função (devido àcrescente oscilação) para os mesmos pontos.

27

4 ConclusõesNo nosso trabalho tirámos as seguintes conclusões acerca do Método das

Soluções Fundamentais (MSF):

• O MSF é um método de convergência muito rápida, a menos de maucondicionamento;

• Os resultados não dependem dos domínios escolhidos;

• A proximidade das fronteiras auxiliares afecta tanto o condicionamentocomo a qualidade da aproximação;

• A introdução da fronteira auxiliar Γ∗2 não melhora os resultados, aocontrário da expectativa teórica;

• O método é muito sensível ao erro;

• O método é sensível à oscilação nas funções dadas;

28

Referências[1] Alves, C.J.S., Chen, C.S.: A new method of fundamental solutions ap-

plied to nonhomogeneous elliptic problems, Advances in ComputationalMathematics (2005) 23 (’05)

[2] Evans, L.C.: Partial Differential Equations, AMS (2005)

[3] Alves, C.J.S.: Apontamentos de Problemas Inversos e Equações Difer-enciais Parciais, IST (2006)

[4] Isakov, Victor: Inverse Problems for Partial Differential Equations,Springer-Verlag New York, Inc. (1998)

[5] Colton, David, Kress, Rainer: Inverse Acoustic and ElectromagneticScattering Theory (Second Edition), Springer-Verlag New York

29

5 Anexo

5.1 ProblemaDirecto.m

30

Documents

André Fernandes Vasconcelos 21 de Julho de 2007calves/AN/lmac/Projecto/... · bem como aperfeiçoar os métodos numéricos que nos permitem prever a na- tureza, localização, e/ou