gusucode.com > demos工具箱matlab源码程序 > demos/brussode.m
function brussode(N) %BRUSSODE Stiff problem modelling a chemical reaction (the Brusselator). % The parameter N >= 2 is used to specify the number of grid points; the % resulting system consists of 2N equations. By default, N is 20. The % problem becomes increasingly stiff and increasingly sparse as N is % increased. The Jacobian for this problem is a sparse constant matrix % (banded with bandwidth 5). % % The property 'JPattern' is used to provide the solver with a sparse % matrix of 1's and 0's showing the locations of nonzeros in the Jacobian % df/dy. By default, the stiff solvers of the ODE Suite generate Jacobians % numerically as full matrices. However, when a sparsity pattern is % provided, the solver uses it to generate the Jacobian numerically as a % sparse matrix. Providing a sparsity pattern can significantly reduce the % number of function evaluations required to generate the Jacobian and can % accelerate integration. For the BRUSSODE problem, only 4 evaluations of % the function are needed to compute the 2N x 2N Jacobian matrix. % % Setting the 'Vectorized' property indicates the function f is % vectorized. % % E. Hairer and G. Wanner, Solving Ordinary Differential Equations II, % Stiff and Differential-Algebraic Problems, Springer-Verlag, Berlin, % 1991, pp. 5-8. % % See also ODE15S, ODE23S, ODE23T, ODE23TB, ODESET, FUNCTION_HANDLE. % Mark W. Reichelt and Lawrence F. Shampine, 8-30-94 % Copyright 1984-2014 The MathWorks, Inc. % Problem parameter, shared with the nested function. if nargin<1 N = 20; end tspan = [0; 10]; y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)]; options = odeset('Vectorized','on','JPattern',jpattern(N)); [t,y] = ode15s(@f,tspan,y0,options); u = y(:,1:2:end); x = (1:N)/(N+1); figure; surf(x,t,u); view(-40,30); xlabel('space'); ylabel('time'); zlabel('solution u'); title(['The Brusselator for N = ' num2str(N)]); % ------------------------------------------------------------------------- % Nested function -- N is provided by the outer function. % function dydt = f(t,y) % Derivative function c = 0.02 * (N+1)^2; dydt = zeros(2*N,size(y,2)); % preallocate dy/dt % Evaluate the 2 components of the function at one edge of the grid % (with edge conditions). i = 1; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + c*(1-2*y(i,:)+y(i+2,:)); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + c*(3-2*y(i+1,:)+y(i+3,:)); % Evaluate the 2 components of the function at all interior grid points. i = 3:2:2*N-3; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ... c*(y(i-2,:)-2*y(i,:)+y(i+2,:)); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ... c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:)); % Evaluate the 2 components of the function at the other edge of the grid % (with edge conditions). i = 2*N-1; dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + c*(y(i-2,:)-2*y(i,:)+1); dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + c*(y(i-1,:)-2*y(i+1,:)+3); end % ------------------------------------------------------------------------- end % brussode % --------------------------------------------------------------------------- % Subfunction -- the sparsity pattern % 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 % ---------------------------------------------------------------------------