PVC linear

Um problema de valor de contorno (PVC) é uma equação diferencial com condições adicionais prescritas não todas em um único ponto. Nesta aula, vamos considerar o problema de valor de contorno de segunda ordem \begin{align} \label{pvc} y'' &= f(x,y,y'), \quad a \le t \le b,\\[2mm] y(a) &= \alpha, \nonumber\\ y(b) &= \beta. \nonumber \end{align} Veja que uma das condições foi dada em $x=a$ e outra foi fornecida em $x=b.$ Condições prescritas sobre o valor da função são ditas condições de Dirichlet. Podemos ter condições também sobre o valor da derivada da função, chamadas condições de Neumann ou mesmo condições mistas.

Vamos buscar aproximações para os valores da função $y$ sobre 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) = f(x_j,y(x_j),y'(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}. $$

Vamos nos concentrar agora no caso mais simples de de equações lineares, ou seja, quando a equação em \eqref{pvc} é dada por \begin{equation}\label{pvclinear} y'' + p(x) y' + q(x) y = r(x), \end{equation} para funções $p,$ $q$ e $r$ conhecidas. Neste caso, as equações \eqref{disc1}, já usando as aproximações de diferenças finitas acima, ficam $$ { y_{j-1}-2y_j+y_{j+1}\over h^2} + p(x_j) {y_{j+1} - y_{j-1}\over 2h} + 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,$ teremos $a_1y_0.$ Como $y_0$ já é conhecido, este termo passa à direita. O mesmo ocorre na última equação, quando $j=n,$ onde haverá 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 essa 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 não nulo. Para $h$ pequeno o suficiente, é possível garantir que não há nem mesmo a necessidade de pivoteamento parcial.

Este método para resolver um problema de valor de contorno é conhecido como método de diferenças finitas .

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 regular. O bloco de código abaixo faz isso para o PVC $$ xy'' + (3+x^2)y'+3x^2y = 0,\; 1\lt x \lt 3,\quad y(1) =1,\quad y(3)=0.2. $$ 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{pvclinear}.

> ### Funções que definem a EDO no formato:
> ### y'' + p(x) y' + q(x) y = r(x)
> p = @(x) (3+x.^2) ./ x;
> 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
> alpha = 1;              # valor de y no extremo esquerdo
> beta = 0.2;             # 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) = alpha;           # 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 motar 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 sua 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. 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.

2. 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. $$