PVC linear

Na aula passada, lidamos com o PVC \begin{align*} y'' &= f(x,y,y'), \quad a \lt x \lt b,\\[2mm] y(a) &= \alpha,\\ y(b) &= \beta. \end{align*} Nesta aula vamos nos concentrar no problema de valor de contorno linear com condições de fronteira dadas sobre os valores de função (conhecidas como condições de Dirichlet), isto é, na situação em que $f(x,y,y') = p(x)y'(x) + q(x)y(x)-r(x)$. Assim, nosso problema alvo é \begin{align} \label{pvc} - y'' + p(x)y' + q(x)y &= r(x), \quad a \lt x \lt b,\\[4pt] y(a)& = \alpha \nonumber \\ y(b)& =\beta. \nonumber \end{align}

A solução deste problema é uma função $y:[a,b]\to\R$. O método numérico que consideraremos nesta aula não tenta aproximar essa função, mas tão somente aproximar amostras dessa função em uma malha regular para o intervalo $[a,b].$ Mais claramente, dado $n\in\mathbb{N},$ seja $h=(b-a)/n$ e defina $x_j = a + jh.$ Queremos encontrar $y_j,$ tal que $$ y_j \approx y(x_j), \quad j=0,1,\ldots,n.$$ Conhecidos os valores $y_j,$ se realmente quisermos uma função contínua que aproxime a função $y,$ podemos resolver um problema de interpolação por partes. Se quisermos uma função diferenciável, podemos realizar a interpolação por splines cúbicas.

Como temos $(n+1)$ incógnitas a determinar, os valor $y_j,$ precisamos de $(n+1)$ equações. As condições de contorno de Dirichlet fornecem duas equações que imediatamente determinam os valores de $y$ nos extremos do intervalo. Com efeito, \begin{align*} y_0 & \equiv \alpha = y(x_0),\\ y_n & \equiv \beta = y(x_n). \end{align*} Nos resta conseguir outras $(n-1)$ equações envolvendo $y_j,$ $j=1,\ldots,(n-1).$ Como a equação diferencial \eqref{pvc} vale para todo $x\in(a,b),$ em particular vale nos pontos $x_j$ da malha, ou seja, \begin{equation}\label{disc1} -y''(x_j) + p(x_j)y'(x_j) + q(x_j)y(x_j) = r(x_j), \quad j=1,\ldots,(n-1). \end{equation} As equações que precisamos são obtidas das equações acima quando $y(x_j)$ é substituído por $y_j$ e $y'(x_j)$ e $y''(x_j)$ são substituídos por aproximações de diferenças finitas. Relembrando, as aproximações de segunda ordem para derivada primeira e segunda são, respectivamente, $$ y'(x_j) \approx {y_{j+1} - y_{j-1}\over 2h}, \quad y''(x_j) \approx { y_{j-1}-2y_j+y_{j+1}\over h^2}. $$

Ao subtituir essas aproximações de diferenças finitas nas equações \eqref{disc1}, ficamos com $$ -\left({ y_{j-1}-2y_j+y_{j+1}\over h^2}\right) + p(x_j) \left({y_{j+1} -y_{j-1}\over 2h}\right) + q(x_j) y_j = r(x_j), \quad j=1,\ldots,(n-1). $$ Agrupando os termos e multiplicando por $h^2,$ ficamos com $$\overbrace{\left(-1-{h\over 2}p_j\right)}^{a_j}y_{j-1} \overbrace{+ \left(2+h^2q_j\right)}^{b_j}y_j +\overbrace{\left(-1+{h\over 2}p_j\right)}^{c_j}y_{j+1} = h^2r_j, \quad j=1,\ldots,(n-1), $$ onde $p_j = p(x_j),$ $q_j=q(x_j)$ e $r_j=r(x_j).$ Repare que na primeira destas equações, quando $j=1,$ temos $a_1y_0.$ Como $y_0$ já é conhecido, este termo passa à direita. O mesmo ocorre na última equação, quando $j=n,$ onde há o termo $c_{n-1}y_n.$

Finalmente chegamos ao sistema linear $Ay = v,$ onde $$A=\left(\begin{array}{ccccc} b_1&c_1&&&\\ a_2&b_2&c_2&&\\ & \ddots&\ddots&\ddots&\\ &&a_{n-2}&b_{n-2}&c_{n-2}\\ &&&a_{n-1}&b_{n-1} \end{array}\right) \qquad v = \left(\begin{array}{c} h^2r_1-a_1y_0\\\ h^2r_2\\\ \vdots\\\ h^2r_{n-2}\\ h^2r_{n-1}-c_{n-1}y_n \end{array}\right) $$ Uma matriz com esta estrutura é dita tridiagonal. A resolução deste sistema pelo processo de eliminação de Gauss é muito eficiente, uma vez que abaixo da diagonal em cada coluna há apenas um elemento possivelmente não nulo. Para $h$ suficientemente pequeno, é possível garantir que não há nem mesmo a necessidade de pivoteamento parcial.

Resolução numérica no Octave

O Octave não tem uma função que resolve imediatamente um PVC. Devemos construir o sistema linear do método de diferenças finitas e resolvê-lo como um sistema linear como qualquer outro. O bloco de código abaixo faz isso para o PVC $$ -y'' +(12-x^2)y'+3xy = 0,\; 1\lt x \lt 3,\quad y(1) = 4,\quad y(3)=5. $$ Primeiro, definimos as funções que descrevem o PVC e os dados para as condições de contorno. Lembre que estamos assumindo que o PVC está no formato \eqref{pvc}.

### Funções que definem a EDO no formato:
### - y'' + p(x) y' + q(x) y = r(x)
p = @(x) (12-x.^2);
q = @(x) 3*x;
r = @(x) 0*x;

### Condições de contorno
a = 1;                  # extremo esquerdo do intervalo de integração
b = 3;                  # extremo direito do intervalo de integração
alfa = 4;               # valor de y no extremo esquerdo
beta = 5;               # valor de y no extremo direito

Passamos à definição da malha de diferenças finitas e do vetor y que abrigará a solução numérica. Repare que o vetor y já é inicializado com os valores da condição de contorno.

n = 20;                 # quantidade de subintervalos
x = linspace(a,b,n+1)'; # pontos da malha
h = x(2) - x(1);        # passo de discretização
y = zeros(n+1,1);       # vetor solução
y(1) = alfa;            # condição de contorno à esquerda
y(end) = beta;          # condição de contorno à direita

Por fim, construímos os vetores que representam as três diagonais da matriz do sistema linear e os usamos para montar a matriz A. Para isso, usamos a função diag que cria uma matriz com uma única diagonal, a partir de um vetor. O vetor independente v é facilmente computado, mas não podemos esquecer de corrigir suas primeira e última componentes. Neste ponto, para relacionar o código do Octave com a teoria acima, lembre que todos os índices no Octave iniciam-se em $1$ e não em $0.$

aa = -1 - h * p(x)/2;   # subdiagonal inferior de A
bb =  2 + h^2 * q(x);   # diagonal principal de A
cc = -1 + h * p(x)/2;   # subdiagonal superior de A

A = diag(aa(3:end-1),-1) + diag(bb(2:end-1)) + diag(cc(2:end-2),1);

v = h^2 * r(x(2:end-1));
v(1)   = v(1)   - aa(2) * y(1);
v(end) = v(end) - cc(end-1) * y(end);

O sistema é resolvido com o operador \ e o resultado é atribuído ao miolo do vetor y, uma vez que a condição de contorno já foi usada.

y(2:end-1) = A\v;
plot(x,y,'o-')

1. No texto desta aula, foram apresentados alguns blocos de código para resolver no Octave um PVC linear. Converta-os em uma função do Octave. Veja uma sugestão de esqueleto para a função.

function [x,y] = pvclinear(p,q,r,a,b,alfa,beta,n)

# Resolução numérica de PVC linear por diferenças finitas
# -y'' + p(x) y' + q(x) y = r(x), a < x < b
#
# y(a) = alfa
# 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

...

endfunction

function [x,y] = pvclinear(p,q,r,a,b,alfa,beta,n)

# Resolução numérica de PVC linear por diferenças finitas
# -y'' + p(x) y' + q(x) y = r(x), a < x < b
#
# y(a) = alfa
# 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

   x = linspace(a,b,n+1)'; # pontos da malha
   h = x(2) - x(1);        # passo de discretização
   y = zeros(n+1,1);       # vetor solução
   y(1) = alfa;            # condição de contorno à esquerda
   y(end) = beta;          # condição de contorno à direita

   aa = -1 - h * p(x)/2;   # subdiagonal inferior de A
   bb =  2 + h^2 * q(x);   # diagonal principal de A
   cc = -1 + h * p(x)/2;   # subdiagonal superior de A

   A = diag(aa(3:end-1),-1) + diag(bb(2:end-1)) + diag(cc(2:end-2),1);

   v = h^2 * r(x(2:end-1));
   v(1)   = v(1)   - aa(2) * y(1);
   v(end) = v(end) - cc(end-1) * y(end);

   y(2:end-1) = A\v;

endfunction 

2. Resolva numericamente o PVC $$ -y'' - x^3y'+\cos(x) y = -x-2,\; 0\lt x \lt 2,\quad y(0) = 0,\quad y(2) = 6. $$ Exiba o gráfico da solução.

3. Considere o problema de valor de contorno $$ y''+2y=-x, \; 0\lt x\lt 1, \qquad y(0)=0, \quad y(1) =0, $$ cuja solução analítica é $$ y(x) = 2 \sum_{n=1}^\infty \frac{(-1)^{n+1}\sin(n\pi x)}{(n^2\pi^2-2)n\pi}. $$

1. Esboce o gráfico de $y,$ para $x\in[0,1].$
2. Estime $$ I = \int_0^1 y(x) dx. $$