%MONTECARLO last updated 3/20/01 % An interactive routine to approximate area under a curve y = f(x) % using a Monte Carlo estimation procedure. % %PROCEDURE: A curve y = f(x) is sketched on an interval [a,b] over which % it is assumed to be nonnegative. A rectangle with base on the % x-axis over [a,b] is constructed with height d >= max(f(x) % over [a,b]. A stream of randomly generated points in the rectangle % is generated and plotted. A count of the number of points with % y-coordinate less thatn the corresponding y-coordinate on f(x) is % tallied. The approximation of the area under f(x) on [a,b] is computed % as (#tallied points/total # used) times (the area of enclosing rectangle). % %RESTRICTIONS: ==> the function must be in terms of variable x; % no check is made % ==> the routine assumes it is >= 0 on interval [a,b]; % no check is made % %INPUTS: The user is asked to enter the following information and click on % "ACCEPT" buttons. % % <><> click in the input areas prior to trying to type % in the corresponding information <><> % % The information can be changed prior to clicking % the "ACCEPT" button. If the information is not acceptable to the % program defaults are used. At any time you can get HELP, RESTART, or % QUIT. % % ==> formula in x for function f(x) (default is f(x) = x) % ==> interval over which f(x) in >= 0 (default is [0,1]) % ==> d = height of enclosing rectangle (default is max(f(x)) on [a,b]) % ==> number of random trials (default is 100) % %OUTPUT: A sketch of y = f(x) in the rectangle from a<= x < = b, 0 <= f(x) <= d % together with random points. The screen is labeled with the formula % for f(x), [a,b], [0,d], number of trials, number of points with y-coordinate % less that the corresponding value of f(x), and the area approximation. % %OPTIONS: The number of trials used can be re-entered to make further approximations. % To change any other information use the RESTART button. % % % % By: David R. Hill, Mathematics Dept., Temple Univ. % Philadelphia, PA. 19122 Email: dhill001@temple.edu %======================================================================= %CONSTRUCTED for use with NSF Project DEMOs with POSITIVE IMPACT %======================================================================= %strings s0=' '; header='Monte Carlo Area Approximation'; cont='Press ENTER to continue.'; mess1='Enter an expression for f(x) ==>'; mess2='Then Press ACCEPT f(x) Button.'; mess3='Enter interval [a,b] for x-values ==>'; mess4='f(x) must be >= 0 on this interval.'; mess5='Then Press ACCEPT [a,b] Button.'; mess6='Enter d >= max(f(x)) over [a,b] ==>'; mess7='Then Press ACCEPT d Button.'; mess8='ERROR: d set to max of f(x).'; mess9='Enter number of random trials ==>'; mess10='Then Press ACCEPT Trials Button.'; mess11='Improper choice; set to 100 Trials.'; mess12='OPTION: new number of trials. ==>'; mess13='For any other changes click RESTART.'; mess14='If done, click QUIT.'; menu=[ 'Options; '; ' '; '1. Get information on this routine. '; ' '; '2. Execute this routine. '; ' '; '0. QUIT! '; ' ']; %INTRO goodch='N'; while goodch=='N' clc,disp(header),disp(s0),disp(menu) ch=input('Make your choice ==> '); if isempty(ch)==1,ch==99;end if ch==1,help montecarlo,disp(s0),disp(cont),pause,goodch='Y';end if ch==2,goodch='Y';end if ch==0,clc,disp('Montecarlo is over!');goodch='Y';return,end if goodch=='N',disp('IMPROPER choice; try again.');disp(cont);pause;end end %seeding the random number generator rand('state',sum(100*clock)) % % %COLOR settings bkgr='white'; %background color bkgray1=[.8 .8 .8]; %background for editable text %dummy callback donothing=' '; % %CALL back for quit button done = ['close(gcf);clc,disp(''MONTECARLO is over!'')']; %callback for the help button helps='set(gcf,''visible'',''off'');clc,help montecarlo,disp(cont),'; helps=[helps 'pause,set(gcf,''visible'',''on'');']; % %callback for funcin % getfunc='delete(funcin),delete(funcinfrm),set(txtHndl,''string'',mess1),'; getfunc=[getfunc 'set(txtHndl2,''string'',mess2);']; getfunc=[getfunc 'fbox=uicontrol(''Style'',''edit'',''units'',''normal'',']; getfunc=[getfunc '''fontsize'',12,''position'',[.55 .2 .3 .05],']; getfunc=[getfunc '''back'',bkgray1,''fontweight'',''bold'',''call'',donothing);']; getfunc=[getfunc 'fenter=uicontrol(''Style'',''push'',''units'',''normal'',']; getfunc=[getfunc '''fontsize'',12,''position'',[.65 .1 .13 .05],']; getfunc=[getfunc '''string'',''ACCEPT f(x)'',']; getfunc=[getfunc '''back'',bkgray1,''fontweight'',''bold'',''call'',gotfunc);']; % %got function callback gotfunc='set(txtHndl,''string'',s0),set(txtHndl2,''string'',s0);'; gotfunc=[gotfunc 'f=get(fbox,''string'');if isempty(f)==1,f=''x'';end;']; gotfunc=[gotfunc 'fname=[''f(x) = '' f];subplot(datasub);']; gotfunc=[gotfunc 'text(0,.9,fname,''color'',''red'',''fontsize'',14,']; gotfunc=[gotfunc '''fontweight'',''bold'',''erasemode'',''none'');']; gotfunc=[gotfunc 'delete(fbox),delete(fenter),']; gotfunc=[gotfunc 'set(txtHndl,''string'',mess3),set(txtHndl2,''string'',mess4);']; gotfunc=[gotfunc 'set(txtHndl3,''string'',mess5);']; gotfunc=[gotfunc 'abbox=uicontrol(''Style'',''edit'',''units'',''normal'',']; gotfunc=[gotfunc '''fontsize'',12,''position'',[.55 .2 .15 .05],']; gotfunc=[gotfunc '''back'',bkgray1,''fontweight'',''bold'',''call'',donothing);']; gotfunc=[gotfunc 'abenter=uicontrol(''Style'',''push'',''units'',''normal'',']; gotfunc=[gotfunc '''fontsize'',12,''position'',[.55 .1 .15 .05],']; gotfunc=[gotfunc '''string'',''ACCEPT [a,b]'',']; gotfunc=[gotfunc '''back'',bkgray1,''fontweight'',''bold'',''call'',gotab);']; %gotab callback gotab='set(txtHndl,''string'',s0),set(txtHndl2,''string'',s0);set(txtHndl3,''string'',s0);'; gotab=[gotab 'xint=get(abbox,''string'');if isempty(xint)==1,xint=''[0,1]'';end;']; gotab=[gotab 'xval=str2num(xint);a=xval(1);b=xval(2);']; gotab=[gotab 'abname=[''on interval '' xint];subplot(datasub);']; gotab=[gotab 'text(0,.8,abname,''color'',''red'',''fontsize'',14,']; gotab=[gotab '''fontweight'',''bold'',''erasemode'',''none'');']; gotab=[gotab 'delete(abbox),delete(abenter),']; gotab=[gotab 'vecf=vectorize(f);x=a:.01:b;y=eval(vecf);xdata=x;ydata=y;']; gotab=[gotab 'pltsub=subplot(''position'',[.05 .3 .60 .60]);']; gotab=[gotab 'plot(xdata,ydata,''-m'',''linewidth'',2);']; gotab=[gotab 'axis([a b 0 max(y)+.1]);axis(axis),hold on;']; gotab=[gotab 'set(txtHndl,''string'',mess6),set(txtHndl2,''string'',mess7);']; gotab=[gotab 'dbox=uicontrol(''Style'',''edit'',''units'',''normal'',']; gotab=[gotab '''fontsize'',12,''position'',[.55 .2 .15 .05],']; gotab=[gotab '''back'',bkgray1,''fontweight'',''bold'',''call'',donothing);']; gotab=[gotab 'denter=uicontrol(''Style'',''push'',''units'',''normal'',']; gotab=[gotab '''fontsize'',12,''position'',[.55 .1 .15 .05],']; gotab=[gotab '''string'',''ACCEPT d'',']; gotab=[gotab '''back'',bkgray1,''fontweight'',''bold'',''call'',gotd);']; %got d callback gotd='dst=get(dbox,''string'');'; gotd=[gotd 'if isempty(dst)==1,dst=''0'';end;']; gotd=[gotd 'dval=str2num(dst);if dval