28.12 O seguinte sistema é um exemplo clássico de EDOs rígidas que podem ocorrer na cinética de reações químicas.
d c1d t & =-0,013 c1-1000 c1 c3 d c2d t & =-2500 c2 c3 d c3d t & =-0,013 c1-1000 c1 c3-2500 c2 c3
Resolva essas reações de t=0 a 50 com condições iniciais c1(0)= c2(0)=1 e c3(0)=0. 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
d c1d t &= -0,013,c1 - 1000,c1,c3, d c2d t &= -2500,c2,c3, d c3d t &= -0,013,c1 - 1000,c1,c3 - 2500,c2,c3,
com condições iniciais c1(0)=1, c2(0)=1, c3(0)=0, no intervalo de tempo tin[0,50].
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 45) e um método rígido (ode 23s), 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.)

Exercícios de Livros Relacionados