Tutorial super rápido do Octave

Em 2004 eu escrevi um tutorial sobre o MuPAD (inclusive, dei um curso sobre ele). O MuPAD era um software de computação algébrica gratuito, com ótimos recursos. O MuPad foi comprado pela Mathworks e virou o pacote simbólico do Matlab.

O tempo passou, eu vim pra Unicamp, que assina o Mathematica, e eu fiz um tutorial (muito rápido) de uso do Mathematica, voltado para alunos de meus cursos, principalmente do curso de equações diferenciais.

O problema é que o Mathematica é pago. Ele funciona muito bem aqui na Unicamp, mas se a pessoa não tiver num lugar que o assina, vai ser complicado. Daí eu fiz um tutorial, bastante detalhado, do uso do Python para trabalhar com matemática. O problema do Python é que ele não é, de fato, um software matemático (apesar de fazer bem o papel).

Então, aqui está um tutorial para o uso básico do Octave. O Octave é um clone gratuito do Matlab, com sintaxe muito parecida. Para fins de matemática “básica”, ele é um ótimo substituto do Matlab. Ele não é tão poderoso quanto o Mathematica para cálculos simbólicos, mas o pacote symbolic até que dá conta de muita coisa. Hoje em dia existe até uma versão online do Octave, que é muito boa.

A ideia é trabalhar no Octave como em uma calculadora matemática. Quase tudo que envolve números, gráficos e aproximações pode ser feito sem o symbolic. O pacote é útil quando queremos respostas exatas, expressões algébricas, derivadas e integrais simbólicas.

Os comandos abaixo foram escolhidos para serem digitados diretamente no console. Pressione Enter para executá-los. Um ponto e vírgula (;) ao final apenas esconde o resultado na tela.

Para uma tarefa mais longa, o melhor é criar um script e executar todo o código de uma vez.

Se quiser fazer cálculos simbólicos, carregue o pacote no início da sessão:

pkg load symbolic

Se o pacote não estiver instalado, use sudo apt install octave-symbolic (no terminal do sistema, não no console do Octave).

1. Contas básicas e funções

Sem nenhum pacote:

1 + 1
2^10
sqrt(2)
pi
sin(pi/2)

format long           % mostra mais casas decimais
sqrt(2)
format short          % volta à exibição usual

Para atribuir valores:

a = 3;
a + 2

Para definir uma função e avaliá-la, a notação @(x) é muito conveniente:

f = @(x) x.^2 + 2*x + 1;
f(2)
f(3)

Nesse exemplo, f(2) resulta em 9. A função f pode ser utilizada novamente em outras contas, sem criar um arquivo.

Importante: talvez a sintaxe que diferencia o Octave/Matlab do Maple/Mathematica seja essa: você precisa escrever .^, .* e ./ para operações elemento a elemento. Isso faz a mesma fórmula funcionar tanto com um número quanto com um vetor.

Por exemplo, defina duas matrizes

A=[1,2; 3,4]
B=[5,0; 2,1]

O comando A*B é o produto usual de matrizes. Porém, se você quiser obter a matriz que é o resultado da multiplicação entrada a entrada de A por B, então pode usar A.*B.

Você pode avaliar a função f numa lista de números, fazendo o comando:

f([0 1 2 3])

Para substituir valores em funções e expressões, é preciso atenção. Se for uma função, a notação é a usual:

f = @(x) sin(x) + x.^2;
f(2)

Se a expressão for simbólica:

pkg load symbolic
syms x
f = sin(x) + x^2;
subs(f, x, 2)               % resultado simbólico
double(subs(f, x, 2))       % número decimal

Substituição de mais de uma variável:

syms x y
g = x^2 + x*y + y^2;
subs(g, [x y], [2 3])

O resultado é 19. Para converter uma expressão simbólica em função numérica reutilizável:

syms x
f = sin(x) + x^2;
h = function_handle(f);
h(2)

2. Raízes de polinômios

Para resolver numericamente

\[ x^2+x-1=0, \]

basta passar a lista dos coeficientes do polinômio, em ordem decrescente de grau:

roots([1 1 -1])

Para um polinômio de grau três:

p = [1 -6 11 -6];
roots(p)

As raízes são 1, 2 e 3 (a ordem de apresentação pode variar). O mesmo vetor de coeficientes serve para avaliar o polinômio:

polyval(p, 4)

Quando quiser a fórmula exata, use o pacote symbolic:

pkg load symbolic
syms x
solve(x^2 + x - 1 == 0, x)

O resultado corresponde a \((-1\pm\sqrt 5)/2\). Para convertê-lo em números decimais:

s = solve(x^2 + x - 1 == 0, x);
double(s)

O comando solve pertence ao pacote symbolic; roots já faz parte do Octave básico.

3. Uma raiz de uma equação não linear

Para encontrar uma raiz real de

\[ e^x-x-2=0, \]

use:

fzero(@(x) exp(x) - x - 2, [0 2])

Não é preciso definir uma função antes: a notação @(x) permite escrever a expressão diretamente no comando. O intervalo [0 2] tem valores de sinais opostos nas extremidades, permitindo localizar uma raiz. Essa equação também tem uma raiz negativa:

fzero(@(x) exp(x) - x - 2, [-3 -1])

Um exemplo envolvendo uma função trigonométrica:

fzero(@(x) x.^2 + cos(x) - 3, [1 2])

O resultado é aproximadamente 1.7951166418. Também é possível fornecer somente um valor inicial, por exemplo fzero(@(x) exp(x)-x-2, 1), mas indicar um intervalo com mudança de sinal costuma ser mais seguro.

Com o symbolic, podemos escrever a equação diretamente, sem @(x):

pkg load symbolic
syms x
solve(x^2 + x == 3, x)

O comando solve procura uma solução simbólica e pode não fornecer uma expressão explícita para equações transcendentes, como x^2 + cos(x) == 3. Nesses casos, fzero é uma boa alternativa numérica. Observe a diferença entre == (equação) e = (atribuição).

4. Sistemas lineares — A\b

Para resolver

\[ \begin{cases} 2x+y=5,\\ x-y=1, \end{cases} \]

escreva simplesmente:

A = [2 1; 1 -1];
b = [5; 1];
A\b

O resultado é o vetor coluna [2; 1]. Essa maneira funciona da mesma forma para sistemas maiores, desde que representados por matrizes.

É preferível usar A\b a calcular explicitamente inv(A)*b.

Se as entradas forem simbólicas, o mesmo operador também pode trabalhar com matrizes simbólicas:

pkg load symbolic
syms a
A = [a 1; 1 -1];
b = [5; 1];
A\b

Outra opção simbólica, com equações escritas diretamente, é:

syms x y
[sx, sy] = solve(2*x + y == 5, x - y == 1, x, y)

5. Matrizes e operações básicas de álgebra linear

Para definir uma matriz, use colchetes: espaços (ou vírgulas) separam as colunas e o ponto e vírgula separa as linhas.

A = [2 1; 1 3]
B = [1 2; 0 1]
v = [3; 4]                % vetor coluna

Operações básicas:

A + B                   % soma
A*B                     % produto matricial
A*v                     % matriz vezes vetor
A^2                     % quadrado da matriz
A'                      % transposta conjugada (para matrizes reais, transposta)
A.'                     % transposta sem conjugação
inv(A)                  % inversa (se existir)
det(A)                  % determinante
trace(A)                % traço
rank(A)                 % posto
rref(A)                 % forma escalonada reduzida
size(A)                 % dimensões da matriz

Para obter autovalores e autovetores:

eig(A)                  % autovalores
[V, D] = eig(A);        % colunas de V: autovetores; diagonal de D: autovalores
A*V - V*D               % matriz aproximadamente nula

Outras matrizes úteis:

eye(3)                  % identidade 3x3
zeros(2,3)              % matriz nula 2x3
ones(2,2)               % matriz 2x2 com todas as entradas iguais a 1
diag([1 2 3])           % matriz diagonal
diag(A)                 % elementos da diagonal de A

Para resolver \(Ax=b\), use A\b, como vimos na seção 4. Isso costuma ser numericamente preferível a inv(A)*b.

Quando as entradas forem simbólicas, podemos usar sym:

pkg load symbolic
syms a
M = [a 1; 1 a];
det(M)
inv(M)

Exemplo: matriz de transição (cadeia de Markov)

Considere um sistema com dois estados, \(E_1\) e \(E_2\). Em cada etapa, a probabilidade de permanecer em \(E_1\) é \(0{,}9\) e a de permanecer em \(E_2\) é \(0{,}8\). A matriz de transição é

\[ P=\begin{pmatrix} 0.9&0.1\\ 0.2&0.8 \end{pmatrix}. \]

Adotaremos vetores de probabilidade em linha: cada linha de \(P\) soma 1, e a distribuição evolui por \(p_{n+1}=p_nP\). No console:

P = [0.9 0.1; 0.2 0.8];
sum(P, 2)               % verifica que as linhas somam 1
p0 = [1 0];             % começamos certamente no estado E1
p0*P                    % distribuição após uma etapa
p0*P^2                  % após duas etapas
p0*P^20                 % após vinte etapas

Observe o uso de P^20: é a potência da matriz, não a potência de cada entrada. A distribuição se aproxima de \((2/3,1/3)\).

A distribuição estacionária \(\pi\) satisfaz \(\pi P=\pi\) e \(\pi_1+\pi_2=1\). Podemos calculá-la diretamente com um sistema linear:

M = P' - eye(2);
M(2,:) = [1 1];
pi_est = (M\[0; 1])'    % distribuição estacionária, em linha
pi_est*P                % confirma a invariância, salvo erros numéricos

Também podemos examinar os autovalores de \(P\):

eig(P)

Um deles é \(1\), como ocorre para toda matriz de transição finita. Neste exemplo, o outro é \(0,7\).

7. Exponencial de matrizes: expm(A) e expm(A*t)

O Octave calcula a exponencial matricial com expm, sem precisar do pacote symbolic. Sua definição é

\[ e^A=I+A+\frac{A^2}{2!}+\frac{A^3}{3!}+\cdots. \]

Por exemplo:

A = [0 -1; 1 0];
expm(A)                 % calcula e^A
t = pi/2;
expm(A*t)               % calcula e^(At), com t numérico

Nesse caso, \(e^{At}\) é uma matriz de rotação. Para \(t=\pi/2\), o resultado é aproximadamente

\[ \begin{pmatrix}0&-1\\1&0\end{pmatrix}. \]

Cuidado: exp(A) aplica a exponencial a cada entrada da matriz; para calcular \(e^A\), o comando correto é expm(A).

Se quisermos a expressão simbólica de \(e^{At}\), com \(t\) indeterminado:

pkg load symbolic
syms t
A = sym([0 -1; 1 0]);
simplify(expm(A*t))

O resultado é

\[ e^{At}=\begin{pmatrix} \cos t&-\sin t\\ \sin t&\cos t \end{pmatrix}. \]

Essa exponencial é a matriz de transição de estados do sistema linear de EDOs \(u'=Au\): a solução com condição inicial \(u(0)=u_0\) é

\[ u(t)=e^{At}u_0. \]

Por exemplo, numericamente, sem symbolic:

A = [0 -1; 1 0];
u0 = [1; 0];
expm(A*2)*u0            % valor da solução no instante t=2

Para visualizar a trajetória toda no plano, basta usar ode45, já apresentado nas seções de EDOs:

[t, U] = ode45(@(t,u) A*u, linspace(0,2*pi,400), u0);
plot(U(:,1), U(:,2), 'linewidth', 2)
axis equal

8. Sistemas não lineares — fsolve

Para resolver numericamente

\[ \begin{cases} x^2+y^2=5,\\ x-y=1, \end{cases} \]

podemos reunir as incógnitas no vetor u = [x; y]:

F = @(u) [u(1)^2 + u(2)^2 - 5; u(1) - u(2) - 1];
fsolve(F, [1.8; 0.9])

A solução encontrada é aproximadamente (2,1). Há outra solução, (-1,-2), que pode ser obtida mudando a aproximação inicial:

fsolve(F, [-1.2; -1.8])

O fsolve procura uma solução próxima do ponto inicial; não encontra automaticamente todas as soluções.

Se precisar resolver simbolicamente, o comando é solve do pacote symbolic:

syms x y
solve([x^2 + y^2 == 5; x-y == 1], [x; y])

9. Gráficos de funções — ezplot, fplot e plot

Para desenhar \(f(x)=x^2\) sem carregar o symbolic nem definir uma função, digite simplesmente:

ezplot('x^2', [-2 2])

Ou, sem indicar o domínio (o Octave escolhe um intervalo padrão):

ezplot('x^2')

Para plotar \(f(x)=e^x\sin x\):

ezplot('exp(x)*sin(x)', [0 9])

Para limitar o eixo vertical e ativar a grade:

ylim([-200 200])
grid on

Dois ou mais gráficos na mesma janela: hold on conserva o desenho atual, e hold off encerra a superposição.

ezplot('sin(x)', [0 2*pi])
hold on
ezplot('cos(x)', [0 2*pi])
hold off
legend('sin(x)', 'cos(x)')
grid on

Também podemos usar fplot com uma função anônima, sem defini-la anteriormente:

fplot(@(x) exp(x).*sin(x), [0 9])

Se o pacote symbolic estiver carregado, fplot aceita até a expressão simbólica diretamente:

syms x
fplot(x^2, [-2 2])

O comando plot, por outro lado, recebe principalmente dados numéricos. Para usá-lo, construímos um vetor de valores de x e calculamos os respectivos y:

x = linspace(-2, 2, 200);
y = x.^3 - x;
plot(x, y)

10. Curvas de nível, curvas implícitas e gráficos 3D (sem Symbolic)

Considere \(g(x,y)=x^2+2y^2\). É possível produzir as curvas de nível diretamente a partir da expressão entre aspas:

ezcontour('x^2 + 2*y^2', [-3 3], 200)
axis equal

O último argumento (200) controla o número de pontos em cada direção da malha. O padrão de ezcontour é de apenas 60 pontos por direção; um número maior deixa as curvas mais suaves, embora demore um pouco mais.

Para desenhar somente a curva \(g(x,y)=4\), escreva-a como \(g(x,y)-4=0\):

ezplot('x^2 + 2*y^2 - 4', [-3 3], 250)
axis equal

Para a superfície em três dimensões:

ezsurf('x^2 + 2*y^2', [-3 3], 100)

Os comandos ezcontour, ezplot e ezsurf fazem parte do Octave básico. Não é necessário declarar syms x y, nem definir g previamente.

Se quiser especificar exatamente os níveis e melhorar o controle sobre o gráfico, pode construir uma malha:

v = linspace(-3, 3, 200);
[X, Y] = meshgrid(v, v);
Z = X.^2 + 2*Y.^2;

contour(X, Y, Z, [1 2 4 6 8])
figure
contour(X, Y, Z, [4 4])      % somente o nível 4
axis equal
figure
surf(X, Y, Z)

11. Integração — numérica ou simbólica

Uma integral definida pode ser calculada numericamente, sem pacote:

integral(@(x) exp(-x.^2), 0, 1)

Também é possível integrar em intervalos infinitos, quando a integral converge e o integrador consegue calculá-la:

integral(@(x) exp(-x.^2), -Inf, Inf)

Para polinômios, é possível integrar algebricamente mesmo sem symbolic, trabalhando com os coeficientes:

p = [1 0 1];                   % x^2 + 1
q = polyint(p);                 % coeficientes da primitiva
polyval(q, 3) - polyval(q, 0)  % integral de 0 a 3: resultado 12

Para obter a expressão exata de uma integral mais geral, precisamos do symbolic:

pkg load symbolic
syms x
f = x^2*exp(x);
int(f, x)                  % integral indefinida
int(f, x, 0, 1)            % integral definida
int(exp(-x^2), x, -inf, inf)

No último caso, a resposta exata é \(\sqrt\pi\).

12. Derivação

Para calcular a derivada como uma expressão, use symbolic:

pkg load symbolic
syms x
f = x^2*exp(x);
diff(f, x)                % primeira derivada
diff(f, x, 2)             % segunda derivada

Para derivadas parciais:

syms x y
g = x^2*y + sin(x*y);
diff(g, x)
diff(g, y)

Sem symbolic, é possível derivar um polinômio representado por seus coeficientes:

p = [1 -3 2];             % x^2 - 3*x + 2
polyder(p)                 % coeficientes de 2*x - 3

Para dados numéricos amostrados, gradient dá uma aproximação da derivada, não uma fórmula simbólica.

13. Substituir valores em expressões

14. Uma EDO — solução numérica com ode45

Considere o PVI

\[ y'=-y+\sin t,\qquad y(0)=1. \]

Podemos resolver e plotar a solução diretamente:

[t, Y] = ode45(@(t,y) -y + sin(t), [0 10], 1);
plot(t, Y)
grid on

Para estimar o valor de \(y(2)\) a partir dos pontos calculados:

interp1(t, Y, 2)

O ode45 escolhe pontos de tempo por adaptação numérica; interp1 interpola o resultado já calculado. Se quiser um controle maior, consulte help odeset.

Se precisar da solução analítica de uma EDO e ela puder ser obtida simbolicamente:

pkg load symbolic
syms y(t)
dsolve(diff(y,t) == -y, y(0) == 1)

Nesse caso, a resposta é \(e^{-t}\).

15. Sistema linear de EDOs

Vamos usar um exemplo semelhante ao do tutorial do Mathematica:

\[ \begin{cases} x'=-x+y,\\ y'=y-2x,\\ x(0)=1,\quad y(0)=2. \end{cases} \]

Não é necessário escrever um programa ou um arquivo de função. Digite:

F = @(t,u) [-u(1) + u(2); u(2) - 2*u(1)];
[t, U] = ode45(F, [0 2*pi], [1; 2]);

Aqui, u(1) representa \(x\) e u(2) representa \(y\). A matriz U tem uma coluna para cada componente da solução: U(:,1) contém os valores de \(x(t)\), e U(:,2) os de \(y(t)\).

Para os gráficos das duas componentes em função do tempo:

plot(t, U(:,1), t, U(:,2))
legend('x(t)', 'y(t)')
xlabel('t')
grid on

Para a curva solução no plano de fase (a trajetória parametrizada \((x(t),y(t))\)):

figure
plot(U(:,1), U(:,2))
axis equal
grid on
xlabel('x')
ylabel('y')

Se a trajetória aparecer um pouco poligonal, você pode pedir ao ode45 que forneça valores em mais instantes, sem mudar o sistema:

[t, U] = ode45(F, linspace(0, 2*pi, 800), [1; 2]);
plot(U(:,1), U(:,2), 'linewidth', 2)
axis equal

O método continua adaptando internamente seus passos para controlar o erro; os 800 instantes servem para amostrar mais densamente a solução que será desenhada.

Nesse exemplo, a solução exata é

\[ x(t)=\cos t+\sin t,\qquad y(t)=2\cos t. \]

Se quiser tentar obter a solução simbolicamente:

pkg load symbolic
syms x(t) y(t)
eqs = [diff(x,t) == -x+y; diff(y,t) == y-2*x];
s = dsolve(eqs, x(0) == 1, y(0) == 2);
s.x
s.y

A resolução analítica é opcional; para obter e desenhar as trajetórias, ode45 já basta.

16. Sistema não linear de EDOs

Considere agora

\[ \begin{cases} x'=-y+xy^2,\\ y'=x-x^3,\\ x(0)=0,\quad y(0)=0.5. \end{cases} \]

A sintaxe numérica é praticamente a mesma:

F = @(t,u) [-u(2) + u(1)*u(2)^2; u(1) - u(1)^3];
[t, U] = ode45(F, [0 3*pi], [0; 0.5]);

Para visualizar x(t) e y(t):

plot(t, U(:,1), t, U(:,2))
legend('x(t)', 'y(t)')
grid on

Para plotar a trajetória no plano:

figure
plot(U(:,1), U(:,2))
axis equal
grid on
xlabel('x')
ylabel('y')

Para desenhar essa trajetória com mais pontos e linhas mais espessas, pode repetir o cálculo assim:

[t, U] = ode45(F, linspace(0, 3*pi, 800), [0; 0.5]);
plot(U(:,1), U(:,2), 'linewidth', 2)
axis equal

O procedimento funciona para outros sistemas de duas equações e, com o mesmo princípio, para sistemas de dimensão maior. O vetor de valores iniciais deve ter tantas entradas quanto o sistema tem incógnitas.

Para uma EDO de segunda ordem, também podemos transformá-la em um sistema de primeira ordem. Por exemplo, \(y''+y=0\), com \(y(0)=1\), \(y'(0)=0\):

F = @(t,u) [u(2); -u(1)];
[t, U] = ode45(F, [0 4*pi], [1; 0]);
plot(t, U(:,1))             % y(t)

17. Tabela rápida: Mathematica → Octave

18. Algumas dicas finais

• Trabalhe no console: Enter executa o comando; não é preciso criar arquivos.
• Prefira os comandos nativos quando só quiser números: roots, fzero, fsolve, A\b, eig, expm, ezplot, integral e ode45.
• Para expressões e gráficos rápidos, não é preciso criar uma função: fzero(@(x) x.^2-3, [1 2]) e ezplot('x^2', [-2 2]) são comandos completos.
• Use symbolic para expressões exatas: syms, solve, diff, int, subs, dsolve.
• = e == têm funções diferentes: = atribui valor; == expressa uma igualdade (muito comum nos cálculos simbólicos).
• .^, .* e ./ são operações elemento a elemento, fundamentais quando a função recebe vetores ou matrizes.
• figure abre outra janela; hold on mantém curvas nos mesmos eixos; axis equal usa a mesma escala nos eixos x e y.
• clear apaga as variáveis da sessão caso queira começar outra conta e reaproveitar nomes.
• Ajuda rápida: help roots, help fzero, help fsolve, help ezplot, help graphics_toolkit, help ode45, help eig, help expm, help dsolve.

Referências

• Manual do GNU Octave
• Documentação do pacote Symbolic