gusucode.com > demos工具箱matlab源码程序 > demos/ihb1dae.m
function ihb1dae %IHB1DAE Stiff differential-algebraic equation (DAE) from a conservation law. % IHB1DAE runs a demo of the solution of a stiff differential-algebraic % equation (DAE) system expressed as a fully implicit system f(t,y,y') = 0. % % The Robertson problem, as it is coded in HB1ODE, is a classic test problem % for codes that solve stiff ODEs. The problem is % % y(1)' = -0.04*y(1) + 1e4*y(2)*y(3) % y(2)' = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2 % y(3)' = 3e7*y(2)^2 % % It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0 % to steady state. % % These differential equations satisfy a linear conservation law that can % be used to reformulate the problem as the DAE % % 0 = y(1)' + 0.04*y(1) - 1e4*y(2)*y(3) % 0 = y(2)' - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2 % 0 = y(1) + y(2) + y(3) - 1 % % This problem is used as an example in the prolog to LSODI [1]. Though % consistent initial conditions are obvious, the guess y(3) = 1e-3 is used % to test initialization. A logarithmic scale is appropriate for plotting % the solution on the long time interval. y(2) is small and its major % change takes place in a relatively short time. Accordingly, the prolog % to LSODI specifies a much smaller absolute error tolerance on this % component. Also, when plotting it with the other components, it is % multiplied by 1e4. The natural output of the code does not show clearly % the behavior of this component, so additional output is specified for % this purpose. % % [1] A.C. Hindmarsh, LSODE and LSODI, two new initial value ordinary % differential equation solvers, SIGNUM Newsletter, 15 (1980), % pp. 10-11. % % See also ODE15I, ODE15S, ODE23T, ODESET, HB1ODE, HB1DAE, FUNCTION_HANDLE. % Mark W. Reichelt and Lawrence F. Shampine % Copyright 1984-2014 The MathWorks, Inc. % Use an inconsistent initial condition to test initialization. y0 = [1; 0; 1e-3]; yp0 = [0; 0; 0]; tspan = [0 4*logspace(-6,6)]; M = [1 0 0 0 1 0 0 0 0]; % Use the LSODI example tolerances and set the value of df/dy' options = odeset('RelTol',1e-4,'AbsTol',[1e-6 1e-10 1e-6], ... 'Jacobian',{[],M}); % For this problem you cannot fix more than two components of the guesses. % Fixing the first two components of y0 leads to the same consistent initial % conditions found by ODE15S in HB1DAE. [y0,yp0] = decic(@f,0,y0,[1 1 0],yp0,[],options); [t,y] = ode15i(@f,tspan,y0,yp0,options); y(:,2) = 1e4*y(:,2); figure; semilogx(t,y); ylabel('1e4 * y(:,2)'); title('Robertson DAE problem with a Conservation Law, solved by ODE15I'); xlabel('This is equivalent to the stiff ODEs coded in HB1ODE.'); % -------------------------------------------------------------------------- function res = f(t,y,yp) res = [ yp(1) + 0.04*y(1) - 1e4*y(2)*y(3) yp(2) - 0.04*y(1) + 1e4*y(2)*y(3) + 3e7*y(2)^2 y(1) + y(2) + y(3) - 1 ];