28.12 O seguinte sistema é um exemplo clássico de EDOs rígidas que podem ocorrer na cinética de reações químicas.
Resolva essas reações de
a 50 com condições iniciais
e
. Se você tiver acesso ao software MATLAB, use ambas as funções padrão (por exemplo, ode45) e rígida (por exemplo, ode23s) para obter suas soluções.
Passo 1
O enunciado apresenta um sistema de EDOs rígidas típico de cinética química, dado por
com condições iniciais
,
,
, no intervalo de tempo
.
Conteúdo: métodos numéricos para EDOs rígidas. Vamos implementar o sistema em MATLAB e integrá-lo com um método não rígido (
Conteúdo: métodos numéricos para EDOs rígidas. Vamos implementar o sistema em MATLAB e integrá-lo com um método não rígido (
ode
) e um método rígido (ode
s
), comparando resultados e tempos de CPU.Passo 2
Defina uma função MATLAB que retorna o vetor
. Por exemplo crie o arquivo kinetics.m com o conteúdo:
function dc \(=\) kinetics(t,c)
dc\(1=\) -\(0.013\)*c(\(1\)) - \(1000\)*c(\(1\))*c(\(3\));
dc\(2=\) -\(2500\)*c(\(2\))*c(\(3\));
dc\(3=\) -\(0.013\)*c(\(1\)) - \(1000\)*c(\(1\))*c(\(3\)) - \(2500\)*c(\(2\))*c(\(3\));
dc \(=\) [dc\(1\); dc\(2\); dc\(3\)];
end
Passo 3
No script principal, resolva com ode
:
tspan \(=\) [\(050\)];
c\(0=\) [\(1\); \(1\); \(0\)];
tic;
[tt\(45\), cc\(45\)] \(=\) ode\(45\)(@kinetics, tspan, c\(0\));
time\(45=\) toc;
Isso usa um método de Runge–Kutta de passo adaptativo não rígido, que pode demandar muitos passos por causa da rigidez.
Passo 4
Em seguida, resolva com ode
s:
tic;
[tt\(23\)s, cc\(23\)s] \(=\) ode\(23\)s(@kinetics, tspan, c\(0\));
time\(23\)s \(=\) toc;
Esse é um método implícito para EDOs rígidas, geralmente exigindo menos passos.
Passo 5
Compare os resultados e o desempenho:
fprintf('ode(45): %d passos, tempo (=) %g s\n', length(tt(45)), time(45));
fprintf('ode(23)s: %d passos, tempo (=) %g s\n', length(tt(23)s), time(23)s);
figure;
plot(tt(45), cc(45), '--', tt(23)s, cc(23)s, '-');
legend('c(1) ode(45)','c(2) ode(45)','c(3) ode(45)', ...
'c(1) ode(23)s','c(2) ode(23)s','c(3) ode(23)s');
xlabel('t'); ylabel('\(c_i\)(t)');
title('Comparação ode(45) vs ode(23)s');
Você verá que as curvas coincidem (mesma solução física) mas ode(23)s usa bem menos passos e é muito mais rápido.
Resposta
Código MATLAB completo e execução típica:
% kinetics.m
function dc = kinetics(t,c)
dc(1) = -0.013*c(1) - 1000*c(1)*c(3);
dc(2) = -2500*c(2)*c(3);
dc(3) = -0.013*c(1) - 1000*c(1)*c(3) - 2500*c(2)*c(3);
dc = [dc(1); dc(2); dc(3)];
end
% main.m
tspan = [0 50];
c0 = [1; 1; 0];
tic;
[tt45, cc45] = ode45(@kinetics, tspan, c0);
time45 = toc;
tic;
[tt23s, cc23s] = ode23s(@kinetics, tspan, c0);
time23s = toc;
fprintf('ode45: %d passos, tempo (=) %g s\n', length(tt45), time45);
fprintf('ode23s: %d passos, tempo (=) %g s\n', length(tt23s), time23s);
figure;
plot(tt45, cc45, '--', tt23s, cc23s, '-');
legend('\(c1_{45}\)','\(c2_{45}\)','\(c3_{45}\', '\(c1_{23s}\)','\(c2_{23s}\)','\(c3_{23s}\');
xlabel('t'); ylabel('\(c_i\)(t)');
title('ode45 vs ode23s - EDO rígida');
Resultados típicos mostram que ambas as soluções coincidem, mas ode23s resolve muito mais eficientemente (muito menos passos e menor tempo de CPU), confirmando que para sistemas rígidos é recomendado usar um integrador implícito.)