function prob_2_32 % % Written by: % -- % John L. Weatherwax 2005-05-07 % % email: wax@alum.mit.edu % % Please send comments and especially bug reports to the % above email address. % %----- % Define the physical constants: global kappa beta a Cp K ga sinth Phi Ph Af G0 kappa = 0.171446272015689e-8; beta = 0.213024626664637e-2; a = 0.108595374561510e+4; Cp = 0.496941623289027e+4; K = 10; ga = 9.80665; sinth = 1; Phi = 1.1e+5; Ph = 797.318; Af = 3.82760; %G0 = 270.9; tD = linspace( 0, 1, 11 ); G0Array = 270.9*( 0.8 + 0.2*exp(-tD) ); options = odeset('Mass',@mass,'MassSingular','no','Events',@events); zU = []; for g0i=G0Array, G0 = g0i; [z,y,ze,ye,ie] = ode45(@odes,[0 5],[795.5; 255.0],options); zU = [zU, ze(end)]; end plot(tD,zU,'-o'); %=================================================================== function dydz = odes(z,y) global kappa beta a Cp K ga sinth Phi Ph Af G0 rho = y(1); T = y(2); dydz = [ (-K*G0*abs(G0/rho) - rho*ga*sinth) (a^2 *Phi*Ph*kappa)/(Cp*Af) ]; function A = mass(z,y) global kappa beta a Cp K ga sinth Phi Ph Af G0 rho = y(1); T = y(2); A = zeros(2); A(1,1) = 1/(rho*kappa) - (G0/rho)^2; A(1,2) = beta/kappa; A(2,1) = -(a^2 *beta*(T + 273.15)*G0)/(Cp*rho^2); A(2,2) = G0/rho; function [value,isterminal,direction] = events(z,y) isterminal = 1; direction = 0; rho = y(1); T = y(2); rhosat = -3.3*(T - 290.0) + 738.0; value = rho - rhosat;