function [ ] = prob_1_15 % % 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. % %----- % Parameters: global a c; a=2; c=1; % Initial conditions: t(1) = 0; x(1) = 1; y(1) = 3; % Final time: tmax = 10.0; % A sample of step sizes: hArray = 10.^[-1:-1:-4]; %hArray = 10.^-4; for hi=hArray, disp( [ 'hi = ', num2str(hi) ] ); t(2:end) = []; x(2:end) = []; y(2:end) = []; [ t, x, y ] = eulerMethod( hi, t, x, y, tmax ); figure; plot( x, y, '-o' ); plotBigG( t, x, y, hi ); end return; %------------------------------------------------------------------------ %------------------------------------------------------------------------ function [ t, x, y ] = eulerMethod( hi, t, x, y, tmax ) global a c; %while( max(t) < tmax ) while( t(end) < tmax ) t(end+1) = t(end) + hi; %x(end+1) = x(end) + hi*a*x(end)*(1-y(end)); %y(end+1) = y(end) - hi*c*x(end)*(1-y(end)); x(end+1) = x(end) + hi*a*x(end)*(1-y(end)); y(end+1) = y(end) - hi*c*y(end)*(1-x(end)); end return; %------------------------------------------------------------------------ %------------------------------------------------------------------------ function [] = plotBigG( t, x, y, hi ) global a c; G = x.^(-c) .* y.^(-a) .* exp( c*x + a*y ); figure; plot( t, G, '-o' ); xlabel( 'time (s)' ); ylabel( 'G (unitless)' ); title( [ 'hi = ', num2str(hi) ] ); return;