Dada a equação de transferência radiativa $$\mu \frac{\partial}{\partial \tau} I(\tau,\mu) + I(\tau,\mu) = \frac{\varpi}{2} \sum^{L}_{l=0}\beta_l P_l(\mu) \int^{1}_{-1}P_l(\mu')I(\tau,\mu')d\mu' + Q(x, \mu),$$ deseja-se, por exemplo, estimar o parâmetro $\varpi$, chamado de albedo, que indica a razão entre radiação incidente e refletida pela superfície do corpo; ou, o parâmetro $L$ que indica o grau de espalhamento da radiação no meio; como também uma fonte externa, dada pelo termo $Q$. Na equação, $\tau \in (0,\tau_f)$ é a variável espacial, sendo $\tau_0$ a espessura óptica do meio; $\mu \in [-1,1]$ é a variável angular, o cosseno do ângulo polar $\theta$ formado entre a direção do feixe de radiação e o eixo $\tau$; e, $I(\tau,\mu)$ é a intensidade radiativa na posição $\tau$, segundo a direção $\mu$.
Neste relatório, são apresentados resultados de alguns testes de reconstrução do grau de espalhamento $L$ e da fonte externa $Q$, utilizando os métodos ADO (ordenadas discretas analítico) [3], Levenberg-Marquardt [1][4] e gradientes conjugados [2][4] implementados em python.
In [1]:
import numpy as np
from ic import ado
from ic import optim
Estrutura da lista de parâmetros para o método ADO: [ $n$, $N$, [ $F_0(\mu)$, $F_{\tau_f}(\mu)$ ], $\tau_f$, $\varpi$, $L$, $Q$ ].
In [2]:
params_0 = [11, 10, [1., 0.], 5., 0.99, 6, 5]
params_1 = [11, 10, [1., 0.], 1., 0.3, 6, 5]
params_2 = [11, 10, [1., 0.], 10., 0.3, 6, 5]
params_3 = [11, 10, [1., 0.], 1., 0.99, 6, 5]
params_4 = [11, 10, [1., 0.], 10., 0.99, 6, 5]
params_5 = [11, 10, [1., 0.], 50., 0.99, 6, 5] # works for CG.
params_6 = [11, 10, [1., 0.], 5., 0.3, 6, 5]
params_7 = [11, 10, [1., 0.], 50., 0.3, 6, 5]
params = params_0
Os dados experimentais utilizados na reconstrução são gerados a partir dos resultados do problema direto e ruídos aleatórios. São gerados n_dataset conjuntos de dados com ruídos diferentes:
In [3]:
n_dataset = 5
dataset = []
ado_result = ado(*params)
for i in range(n_dataset):
noise = np.random.uniform(0.99, 1.01, params[0])
dataset.append(ado_result * noise)
Estimativas iniciais:
In [4]:
x0 = [3, 0.]
A função a ser minimizada no processo de otimização:
In [5]:
def f(x0):
params[5] = x0[0]
params[6] = x0[1]
return ado(*params)
As estimativas são realizadas utilizando os métodos de Levenberg-Marquardt e dos gradientes conjugados com cada um dos n_dataset conjuntos de dados:
In [6]:
import warnings
tol = 1e-13
max_iter = 70
est_lm = []
est_cg = []
with warnings.catch_warnings():
warnings.simplefilter('ignore')
for data in dataset:
est_lm.append(optim.lm(f, data, x0, tol, max_iter))
est_cg.append(optim.cg(f, data, x0, tol, max_iter))
In [7]:
for i in range(n_dataset):
lm, cg = est_lm[i], est_cg[i]
print('Test %d' % i)
print('------')
print(' L | Q | Iterations')
print('Levenberg-Marquardt: %.10f | %.10f | %d' % (lm[-1][0], lm[-1][1], len(lm)))
print('Conjugate Gradient: %.10f | %.10f | %d\n' % (cg[-1][0], cg[-1][1], len(cg)))
[1] A. J. Silva Neto e F. D. Moura Neto, Problema Inverso: Conceitos Fundamentais e Aplicações, EdUERJ, 2005. v. 1. 82p.
[2] J. R. Shewchuk, An Introduction to the Conjugate Gradient Method Without the Agonizing Pain, School of Computer Science, Carnegie Mellon University, Pittsburgh, PA 15213, August 4, 1994.
[3] L. B. Barichello, Explicit Formulations for Radiative Transfer Problems, In: Helcio R. B. Orlande; O. Fudyin; D. Maillet; R. M. Cotta. (Org.). Thermal Measurements and Inverse Techniques. Boca Raton: CRC Press, 2011, p. 541-562.
[4] M. N. Ozisik e Helcio R. B. Orlande, Inverse Heat Transfer: Fundamentals and Applications, Taylor & Francis, Inc.