Download<< Back
function run_bruss
% Space discretisation (ignoring boundary)
N=100;
x = (1:N)/(N+1);
% Initial conditions:
y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)];
% Time space:
t0 = 0;
T = 100;
% Compute solution
% We parse a sparse matrix (jpattern(N)) which denotes where the
% stiffness matrix is non-zero - this allows ode15s to optimise
% to reduce the number of function evaluations
options = odeset('Vectorized','on','JPattern',jpattern(N));
[t,y] = ode15s(@(t,x) bruss(t,x,N), [t0,T], y0, options);
% Add boundary and split into u and v components
uBC = ones(size(t));
vBC = 3*ones(size(t));
u = [uBC y(:,1:2:end) uBC];
v = [vBC y(:,2:2:end) vBC];
x = [0, x, 1];
% Plot 1D slices in space
mid = round(N/2);
figure;
plot(t, u(:,2), 'k', t, u(:,mid), 'b', t, u(:,N+1), 'r');
legend(['x = ' num2str(x(1))], ['x = ' num2str(x(mid))], ['x = ' num2str(x(N+1))]);
xlabel('t');
ylabel('u');
figure;
plot(t, v(:,2), 'k', t, v(:,mid), 'b', t, v(:,N+1), 'r');
legend(['x = ' num2str(x(1))], ['x = ' num2str(x(mid))], ['x = ' num2str(x(N+1))]);
xlabel('t');
ylabel('v');
% Plot surfaces
figure;
surf(x, t, u);
xlabel('x');
ylabel('t');
zlabel('u');
figure;
surf(x, t, v);
xlabel('x');
ylabel('t');
zlabel('v');
function S = jpattern(N)
% Jacobian sparsity pattern
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end
end
