Download<< Back
function bruss_stiff_estimate
N = 100;
y0 = [1+sin((2*pi/(N+1))*(1:N)); 3*ones(1,N)];
disp('Jacobian')
A = bruss_jacobian(y0,N);
spectral_points = eig(full(A));
% Rescale h = pi/(N+1)
h = pi/(N+1);
spectral_points = spectral_points/h;
% Plot spectrum
figure;
plot(spectral_points,'r*')
title('Spectrum');
xlabel('Re');
ylabel('Im');
% Select (and sort) eigenvalues where the real part > -50
filtered = spectral_points(real(spectral_points)>-50);
[~,index]=sort(real(filtered));
filtered=filtered(index);
figure
plot(filtered,'r*')
title('Rightmost Spectrum');
xlabel('Re');
ylabel('Im');
% Sparsity pattern
figure;
spy(A); % nonzero elements of the stiffness matrix
title('Jacobian sparsity pattern');
zoom on;
function J = bruss_jacobian(y,N)
%BRUSS_JACOBIAN Computes a (sparse) Jacobian at the specified point for Brusselator
c = 0.02 * (N+1)^2;
B = zeros(2*N,5);
B(1:2*(N-1),1) = B(1:2*(N-1),1) + c;
i = 1:2:2*N-1;
B(i,2) = 3 - 2*y(i).*y(i+1);
B(i,3) = 2*y(i).*y(i+1) - 4 - 2*c;
B(i+1,3) = -y(i).^2 - 2*c;
B(i+1,4) = y(i).^2;
B(3:2*N,5) = B(3:2*N,5) + c;
J = spdiags(B,-2:2,2*N,2*N); % Note this is a SPARSE Jacobian.
end
end
