Solução numérica de um PVI no Octave

O problema de valor inicial padrão para o qual vamos construir métodos numéricos é \begin{equation}\label{standard} y'(t) = f(t,y), \quad t>t_0, \qquad y(t_0)=y_0. \end{equation} A função $y$ pode ser vetorial, o que leva $f$ a ser vetorial também, ou seja, $$ y(t) = \begin{pmatrix}y_1(t)\\\ \vdots \\\ y_m(t)\end{pmatrix}, \qquad f(t,y) = \begin{pmatrix}f_1(t,y)\\\ \vdots \\\ f_m(t,y)\end{pmatrix}. $$

Um exemplo de um problema que se encaixa neste formato é \begin{equation}\label{exemplo} y'(t) \equiv \left(\begin{array}{c}y'_1(t)\\[1mm] y'_2(t)\end{array}\right) =\left(\begin{array}{c}ty_2 - 3y_1\\[1mm] -(1-t)^2y_2\end{array}\right) \equiv f(t,y), \quad t\gt 1, \qquad y(1) = \left(\begin{array}{r} 2 \\[1mm] 1\end{array}\right). \end{equation}

Antes de entrar no estudo desses métodos para o problema \eqref{standard}, é preciso definir que tipo de aproximação queremos. Veja que a solução de um problema de valor inicial (PVI) é uma função, logo poderíamos pensar que um método numérico deve fornecer uma função que aproxime a função solução. Há métodos que fazem isso, como os métodos de colocação, por exemplo. Neste curso, entretanto, teremos um objetivo mais simples. Vamos nos contentar em aproximar amostras da função solução, sobre uma malha de pontos.

Considere a malha $t_0\lt t_1 \lt \cdots \lt t_N.$ Estamos interessados em determinar $y_n,$ uma aproximação para $y(t_n),$ o valor da solução de \eqref{standard} em $t_n.$ Se as aproximações $y_0, y_1, \ldots, y_N,$ forem determinadas, definir uma função contínua a partir dessa amostras recai em um problema de interpolação.

No Octave, problemas de valor inicial no formato \eqref{standard} podem ser facilmente resolvidos com a função lsode. Por exemplo, o problema \eqref{exemplo}, seria resolvido com

# Define a função da equação diferencial e
# a condição inicial
f = @(y,t) [t*y(2)-3*y(1); -(1-t)^2*y(2)];
y0 = [2; 1];

# Malha regular em t, onde aproximaremos a solução
t = linspace(1, 3, 21)';

y = lsode(f, y0, t);

plot(t, y(:,1), '.-', t, y(:,2), '.-')
legend("y_1", "y_2")

Repare que na definição da função $f$ no Octave, como exigido pela rotina lsode, a ordem dos argumentos é diferente da que convencionamos aqui. Como a equação diferencial definida pela função $f$ é vetorial, a solução $y$ produzida é uma matriz com duas colunas, sendo a primeira coluna, y(:,1), as aproximações para a $y_1$ nos pontos da malha, enquanto que a segunda coluna, y(:,2), seriam as amostras para a função $y_2.$

Outro exemplo seria o sistema de equações diferenciais de primeira ordem a seguir. \begin{align*} S'(t) & = -\alpha S(t)I(t) + \gamma R(t),\\ I'(t) & = \hphantom{-}\alpha S(t)I(t) - \beta I(t),\\ R'(t) & = \hphantom{-}\beta I(t)-\gamma R(t). \end{align*} Este sistema, conhecido como modelo SIR, é utilizado é um modelo clássico de dinâmica populacional, utilizado na modelagem da evolução de epidemias em uma população . $S(t)$ representa a parcela da população suscetível ao contágio no instante $t,$ a função $I(t)$ representa a parcela da população que está infectada no instante $t$ e $R(t)$ representa a fração da população que está recuperada no instante $t.$ Os parâmetros $\alpha,$ $\beta$ e $\gamma$ quantificam a taxa com que indivíduos em um determinado grupo se movem para outros grupos. Para resolver este problema, basta definir $y(t) = (S(t),I(t),R(t))^T$ e desta forma $f(t,y)$ é dada por $$ f(t,y) = \begin{pmatrix} -\alpha y_1y_2 + \gamma y_3,\\ \hphantom{-}\alpha y_1y_2 - \beta y_2,\\ \hphantom{-}\beta y_2-\gamma y_3. \end{pmatrix}. $$ Por exemplo, a escolha $\alpha=1.7,$ $\beta=0.5$ e $\gamma=0.1,$ com condição inicial $S(0)=y_1(0)=0.95$, $I(0)=y_2(0) = 0.05$ e $R(0)=y_3(0)=0,$ produz o seguinte comportamento.

Por fim, nem todos os problemas de valor inicial estão no formato padrão \eqref{standard}. Na próxima aula veremos como converter um PVI para o formato padrão. Só então começaremos a estudar os métodos numéricos para resolver o problema \eqref{standard} e assim entender o que o comando lsode do Octave pode estar fazendo.

1. Usando o Octave (ou o Octave online), resolva o PVI a seguir, aproximando sua solução em uma malha regular com 31 pontos no intervalo $[0,3].$ Qual o valor estimado para $y(1.5)$? E para $y(3)$? $$ y'(t) = t^2-y,\quad t>0,\qquad y(0)=1. $$

2. Usando o Octave (ou o Octave online), tente reproduzir o gráfico do modelo SIR, exibido no texto.