www.gusucode.com > demos工具箱matlab源码程序 > demos/ddex2.m
function ddex2 %DDEX2 Example 2 for DDE23. % This example solves a cardiovascular model due to J.T. Ottesen, Modelling % of the Baroflex-Feedback Mechanism With Time-Delay, J. Math. Biol., 36 % (1997) when the peripheral pressure R is reduced exponentially from its % value of 1.05 to 0.84 beginning at t = 600. % % The example shows how the 'Jumps' option is used to inform the solver % about discontinuities in low-order derivatives at points known in advance % (t = 600). Instead of using 'Jumps', this problem can be solved by % breaking it into two pieces: % sol = dde23(@ddex2de,tau,history,[0, 600]); % sol = dde23(@ddex2de,tau,sol,[600, 1000]); % The solution structure SOL on [0, 600] serves as history for restarting % the integration at t = 600. In the second call, DDE23 extends SOL so that % the solution is available on all of [0 1000]. When discontinuities occur % in low-order derivatives at points known in advance, it is better to use % the 'Jumps' option. When discontinuities must be located with event % functions, it is necessary to restart at the discontinuities. % % See also DDE23, DDESET, FUNCTION_HANDLE. % Jacek Kierzenka, Lawrence F. Shampine and Skip Thompson % Copyright 1984-2014 The MathWorks, Inc. % Problem parameters, shared with the nested function p.ca = 1.55; p.cv = 519; p.R = 1.05; p.r = 0.068; p.Vstr = 67.9; p.alpha0 = 93; p.alphas = 93; p.alphap = 93; p.alphaH = 0.84; p.beta0 = 7; p.betas = 7; p.betap = 7; p.betaH = 1.17; p.gammaH = 0; P0 = 93; Paval = P0; Pvval = (1 / (1 + p.R/p.r)) * P0; Hval = (1 / (p.R * p.Vstr)) * (1 / (1 + p.r/p.R)) * P0; history = [Paval; Pvval; Hval]; tau = 4; options = ddeset('Jumps',600); sol = dde23(@ddex2de,tau,history,[0, 1000],options); figure plot(sol.x,sol.y(3,:)) title('Heart Rate for Baroflex-Feedback Mechanism.') xlabel('time t') ylabel('H(t)') % ----------------------------------------------------------------------- % Nested function -- problem parameters provided by the outer function % function dydt = ddex2de(t,y,Z) % Differential equations function for DDEX2. if t <= 600 p.R = 1.05; else p.R = 0.21 * exp(600-t) + 0.84; end ylag = Z(:,1); Patau = ylag(1); Paoft = y(1); Pvoft = y(2); Hoft = y(3); dPadt = - (1 / (p.ca * p.R)) * Paoft + (1/(p.ca * p.R)) * Pvoft ... + (1/p.ca) * p.Vstr * Hoft; dPvdt = (1 / (p.cv * p.R)) * Paoft ... - ( 1 / (p.cv * p.R) + 1 / (p.cv * p.r) ) * Pvoft; Ts = 1 / ( 1 + (Patau / p.alphas)^p.betas ); Tp = 1 / ( 1 + (p.alphap / Paoft)^p.betap ); dHdt = (p.alphaH * Ts) / (1 + p.gammaH * Tp) - p.betaH * Tp; dydt = [ dPadt dPvdt dHdt ]; end % ----------------------------------------------------------------------- end % ddex2