Método de shooting na prática

Nesta aula vamos nos concentrar no método de shooting aplicado aa PVC linear \begin{equation}\label{shooting:pvclinear} Ly \equiv -y''+p(x)y'+q(x)y = r(x), \quad a\lt x \lt b, \qquad y(a)=\alpha, \quad y(b) = \beta. \end{equation}

Vimos na aula passada que a solução deste PVC é $y=u(x) + sv(x)$, com \begin{equation}\label{shooting:s} s =\frac{\beta - u(b)}{v(b)}, \end{equation} onde $u$ e $v$ são as soluções dos PVI's \begin{align} Lu =& r(x), \quad a\lt x \lt b, & u(a)=\alpha, \quad u'(a) = 0,\label{pvi1} \\[6pt] Lv =& 0, \quad a\lt x \lt b, & v(a)=0, \quad v'(a) = 1.\label{pvi2} \end{align}

PVI's no formato padrão

Para aplicar métodos numéricos na resolução dos PVI's \eqref{pvi1} e \eqref{pvi2}, é preciso primeiro que eles estejam no formato padrão de PVI, para o qual os métodos numéricos são criados.

Em uma aula anterior, vimos como converter um PVI de segunda ordem em um sistema de equações diferenciais ordinárias de primeira ordem com condições iniciais. O PVI \eqref{pvi1}, no formato padrão, se escreve como \begin{align}\label{pvi1padrao} U(x) &\equiv \begin{bmatrix}u(x) \\\ u'(x) \end{bmatrix}, & U'(x) &= \begin{bmatrix} U_2 \\\ p(x)U_2+q(x)U_1-r(x) \end{bmatrix},& U(a) &= \begin{bmatrix}\alpha \\\ 0 \end{bmatrix}. \end{align} De maneira análoga, o PVI \eqref{pvi2}, no formato padrão, fica \begin{align}\label{pvi2padrao} V(x) &\equiv \begin{bmatrix}v(x) \\\ v'(x) \end{bmatrix}, & V'(x) &= \begin{bmatrix} V_2 \\\ p(x)V_2+q(x)V_1 \end{bmatrix},& V(a) &= \begin{bmatrix}0 \\\ 1 \end{bmatrix}. \end{align}

Exemplo trabalhado

Para efeito de exemplo, considere o PVC \begin{equation}\label{exemplo} y'' - {y\over 4} = 2x, \quad 0\lt x\lt 2, \qquad y(0) = 1, \quad y(2) = 4. \end{equation} O primeiro PVI associado a esse PVC fica \begin{align}\label{pviu} U(x) &\equiv \begin{bmatrix}u(x) \\\ u'(x) \end{bmatrix}, & U'(x) &= \begin{bmatrix} U_2 \\\ {1\over 4}U_1 -2x\end{bmatrix},& U(0) &= \begin{bmatrix}1 \\\ 0 \end{bmatrix}. \end{align} O segundo PVI associado ao PVC é \begin{align}\label{pviv} V(x) &\equiv \begin{bmatrix}v(x) \\\ v'(x) \end{bmatrix}, & V'(x) &= \begin{bmatrix} V_2 \\\ {1\over 4}V_1 \end{bmatrix},& V(0) &= \begin{bmatrix}0 \\\ 1 \end{bmatrix}. \end{align}

No Octave, a resolução numérica do PVC pelo método de shooting pode ser implementada com utilizando o comando lsode, para a resolução dos PVI's (veja nesta aula como).

a = 0; b = 2;                      # Intervalo [a,b]
alpha = 1; beta = 4;               # Valores de contorno para y
x = linspace(a,b,41)';             # Discretização de [a,b] em 41 pontos

Fu = @(U,x) [U(2); U(1)/4 + 2*x];  # Definição do PVI para u
U0 = [alpha;0];                    # Condição inicial para u

Fv = @(V,x) [V(2); V(1)/4];        # Definição do PVI para v
V0 = [0;1];                        # Condição inicial para v

# Resolução numérica dos dois PVI's
U = lsode(Fu, U0, x);
V = lsode(Fv, V0, x);

s = (beta - U(end,1)) / V(end,1);  # Determinação de s

y = U(:,1) + s * V(:,1);           # Solução do PVC
plot(x,U(:,1),x,V(:,1),x,y)        # Gráfico das soluções dos PVI's e do PVC
legend("u","v","y")

A figura abaixo exibe os gráficos das aproximações obtidas para as funções $u,$ em vermelho, $v,$ em azul, e $y,$ em preto.

A solução exata deste PVC é $$ y(x) = e^{x\over 2}+ \left({4+8b-e^{b\over 2}\over \sinh(b/2)}\right)\sinh(x/2) - 8x, $$ onde $b=2$ é o extremos direito do intervalo de integração. Usando-a, podemos calcular o erro máximo entre a aproximação computada e o valor exato.

# Solução exata do PVC em termos de b, extremo direito do intervalo
yexato = @(x,b) exp(x/2) + (4+8*b-exp(b/2))/sinh(b/2) * sinh(x/2) - 8*x;

norm(y-yexato(x,b),inf)
ans = 8.3262e-08

Observe que o erro máximo, computado através da norma infinito da diferença entre a aproximação numérica e o valor exato, ficou um pouco abaixo de $10^{-7}.$ Isto é consistente com a precisão utilizada na resolução dos PVI's, uma vez que a função lsode trabalha com tolerâncias da ordem de $10^{-8},$ por padrão.

Suponha agora que o extremo direito do intervalo para PVC \eqref{exemplo} estivesse em $b=4$ ao invés de $2.$ Repare que os PVI's em \eqref{pviu} e \eqref{pviv} não seriam alterados. Apenas $s$ depende de $b$.

Para simular este caso, rode novamente o primeiro bloco de código, alterando apenas o valor de $b$ na primeira linha. A figura a seguir, análoga à anterior, exibe como ficam os gráficos de $u,$ $v$ e $y$.

Observe que o valor de $u(b)$ cresce bem rápido. Com consequência, a expressão \eqref{shooting:s} pode sofrer perda dedígitos significativos pela grande diferença em ordem de magnitude entre o valor de $\beta$ e o valor de $u(b).$

De modo geral, se $u(b)$ crescer à medida que $b$ cresce, podemos ter que $$ \mathop{\mbox{fl}}(s) = \mathop{\mbox{fl}}\left({\beta - u(b)\over u(b)}\right) \approx - {u(b)\over v(b)}. $$ Além disso, o cálculo de $y(x)$ também sofre com cancelamento de dígitos significativos, quando $x\approx b$, uma vez que $$ y(x) = u(x) + s v(x) \approx u(x) - {u(b)\over v(b)} v(x) = u(x) - \left({v(x) \over v(b)}\right) u(b). $$

Este fenômeno de perda de acurácia próximo do extremo direito do intervalo de integração pode ser observado no PVC $$ -y'' + 4y = 0, \quad 1\lt x\lt 10, \qquad y(1) = 1, \quad y(10) = 1. $$ A figura abaixo exibe os gráficos da aproximação numérica obtida com o método de shooting, em vermelho, e da solução exata, em preta.

Existem alternativas que tentam mitigar essas dificuldades do método de shooting como, por exemplo, repartir o intervalo de integração $[a,b]$ em diversos subintervalos de modo que, em cada um deles, a função $u$ o crescimento da função $u$ fique controlado. Porém, ao invés de discutir essas alternativas, a partir da próxima aula vamos explorar um método mais robusto e muito mais utilizado para resolver problemas de valor de contorno.

Referências

Hebert B. Keller. Numerical Methods for Two-Point Boundary-Value Problems. Dover, 1992.

1. A partir do bloco de código exibido no texto, crie uma função para o Octave que implemente o método de shooting para um PVI linear, com condições de contorno de Dirichlet. Se precisar, consulte esta aula sobre funções no Octave. Sua função deve ter a seguinte estrutura:

function [x,y,u,v,s] = shooting(p,q,r,a,b,alpha,beta,n=21)

# Método de shooting para aproximar a solução do PVC
# -y''+ p(x)y' + q(x)y = r(x),  a < x < b
#
# y(a) = alpha
# y(b) = beta
#
# Esta função retorna:
# x - vetor regularmente amostrado em [a,b] com n amostras
# y - aproximação para a solução do PVC em x
# u e v - aproximações para as soluções dos PVI's relacionados
# s - escalar tal que y = u + s v

...

endfunction

2. Resolva numericamente os PVI's associados a cada PVC método de shooting e observe o comportamento das funções $u$ e $v$. Para quais desses PVC's podemos ter problemas de acurácia?

1. $-y''-{1\over x}y'+{1\over x^2}y=1,$ $1\lt x \lt 4,$ $y(1)=3,$ $y(4)=0.$
2. $-y''+xy'+(4-x^2)y=0$, $0\lt x\lt 2$, $y(0)=1$, $y(2)=2.$
3. $-y''+xy'+(x-5)^2y=0$, $0\lt x\lt 8$, $y(0)=1$, $y(8)=2.$
4. $-y''-{2\over 1+x^2}y' + {1\over(4+x^2)}y=x,$ $1\lt x\lt b,$ $y(1)=0,$ $y(b)=1,$ com $b = 4,$ $8$ e $16$.
5. $-y''+xy-4y=0,$ $0\lt x\lt b,$ $y(0)=y(b)=1,$ com $b = 3,$ $4$ e $5.$