{VERSION 4 0 "IBM INTEL NT" "4.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Century Schoolb ook" 1 12 0 0 0 0 2 2 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Text Output" -1 2 1 {CSTYLE "" -1 -1 "Courier" 1 10 0 0 255 1 0 0 0 0 0 1 3 0 0 1 }1 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "R 3 Font 0" -1 256 1 {CSTYLE "" -1 -1 "Lucida Sans Typewriter" 1 12 0 0 0 0 2 2 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "R3 \+ Font 2" -1 257 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 0 0 0 2 2 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }} {SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "restart:with(plots):with(stats):read `session5.txt`: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 " UNIT 5 THE RANDOM WALK" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 160 "(HOW TO USE \"ITER\" FOR EULER TYPE PROBLEMS AND THE RANDOM WALK WHER E INPUT CAN BE READ FROM AN ARRAY OR A RANDOM GENERATOR AND OUPUT CAN \+ BE WRITTEN TO AN ARRAY)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 30 "[1] A random number generator:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 659 "In this unit we shall us e an explicit numerical technique to solve an ordinary differential eq uation with a right hand side that is given numerically, rather than a s an analytic function. In fact in the example we choose the right han d side of the equation will be a random variable so that it is impossi ble to give it an analytic form. In this case we cannot use dsolve( ,n umeric) and must write the numerical scheme ourselves. We begin by lea rning how to generate random numbers in Maple. Entering the following \+ command will make available a generator of numbers distributed normall y with mean zero and standard deviation 1 , given to 2 significant fig ures." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 2 "" 0 "" {TEXT -1 1 " \n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "N:=stats[random,normald[0,1]] (generator[3]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 64 "To get a random number from this s equence use 'N()'$1; Try this:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "N()$1;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 90 "The next command can be us ed to display any number of values drawn from this distribution." }} {PARA 0 "" 0 "" {TEXT -1 96 "Edit the command to obtain some further s equences of random numbers from an N(0,1) distribution." }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 17 "seq(N(),i=1..10);" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 425 "We would l ike to check that this is a normal distribution. To do this we count t he number of times the output falls in bins of width .01 in a long seq uence of trials. This is achieved in the second line of code below. Wi th a little thought you should be able to explain how this line does t he counting.(The first line initialises the counters and the third lin e sets up a list for plotting: these bits need not concern you.) " }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 181 "l:=[]:for k from 1 to 100 do n[k]:=0 od:'N()'$1: \nf or k from 1 to 1500 do j:=min(max(1,round((N()*10+50))),100):n[j]:=n[j ]+1 od:\nfor kk from 1 to 100 do l:=[op(l),[kk-50,n[kk]]] od:" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 88 "If you want to inspect the data (not to be recommended if you have a lot of it!) type l;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "The next line plots the data. " }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 30 "plot(l,x=-50..50, style=line);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 116 "Exercise[1] Expla in how the random numbers are being binned in the above Maple code, an d interpret the plot results." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "[2] The Euler Method revisited." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 609 "To rec all how to use ITER to implement the Euler (forward difference) method edit the following commands to solve a first order equation with a ri ght hand side of your own choice (ie of the form y'=f(x,y) where in th e example f(x,y)=-xy). There is no need to make it too complicated! If you choose an example for which Maple cannot obtain an analytical sol ution then ITER will fail. To prevent this use either the ITER(......, noplot) option, which will supress the attempt to plot an analytical s olution or the ITER(.......,numeric) option which will attempt to comp are your numerical method with Maple's own." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "de2:=D(y)(x)+x*y(x)=0; Dy:=(x,y)->-x*y; y(a+h): =(x,y,stepsiz)->y+stepsiz*Dy(x,y); " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "x(a+h):=(x,stepsiz)->x+stepsiz;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "ITER(de2,x(a+h),y(a+h),npoints=20,step=.05,xs tart=0,ystart=1,noplot);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 " " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 146 "Exercise [2] Use the ITER(.....,numeric) option to \+ compare your Euler method with Maple's numerical integration routine a nd comment on the result." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "[3] Brownian moti on: A stochastic differential equation" }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}{PARA 0 "" 0 "" {TEXT -1 255 "We shall now use the Euler method to solve a differential equation for the motion of a particle under a ra ndom force, y'+y=F(x), where F is drawn from a normal distribution N(0 ,1). First we set up the normal distribution rounded to 2 significant \+ figures. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "N:=stats[random,normal d[0,1]](generator[3]):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 44 "Then we define the procedure to obtain F(x)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "F:=(x)->N();" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 235 "Exercise [3] Now copy and edit the text from section [2] to obtain the solution to dy/dx = F(x)-y using ITER. Start from y=0 a t x=0. Use a reasonable number of points! Do not let Maple attempt to \+ find a comparison solution (it can't!)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 561 "The equation you have just solved describes the velocity of a \+ particle in Brownian motion (The general form is y'+ay =f (x,y) - wit h x the time and y(x) the velocity where a is a constant (=1 above). I n general f may depend on y and contain a deterministic (non-random) c omponent. This general form is called a Langevin equation - but you do n't need to know that!) The ay term represents the effect of collision s in reducing the momentum of the particle (viscosity) and f (the driv ing force) represents the effect of collisions inducing the Brownian \+ movement." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 209 "Exercise[4] Investigate the effect of changing the magnitude of t he viscosity, the driving force, and the initial velocity. Paste in an y interesting results (2 to 3 will be quite enough) and comment on th em. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 336 "We may wish to store th e solution at each of the calculated points. ITER has a dataout= opti on to do this. dataout= must be set to a predefined array containing t he appropriate number of points as in the following. Edit the command \+ to work for your number of points and your names for the iteration pro cedures in x (iscx) and y (iscy)!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 89 "veco:=array(1..251);ITER(de,iscx,iscy,npoints=250,step=.1,x1=0,y1= 0,dataout=veco,noplot):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 309 "This is \+ useful in the Brownian motion case because it allows us to integrate t he velocity to obtain the position of the Brownian particle. ITER has \+ a datain= option to do this. The following takes the differential equa tion y'=Z, where Z is the output vector containing the particle veloci ties at each time ,x," }}{PARA 0 "" 0 "" {TEXT -1 160 "and and integra tes it to find the displacements y(x). Note the names of the function s have been changed for convenience to distinguish them from those ab ove. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 297 "Exercise [5] Edit the text to integrate your Brownian motion equation . You may wish to do some further editing to tidy up the names of the \+ procedures and the choice of variables! Paste in the resulting displac ement [together with the corresponding velocity if you have not alread y done so above.] " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "de1:=D(y)(x)= Z; Dy1:=(x,y)->Z; iscy1:=(x,y,stepsiz)->y+stepsiz*Dy1(x,y); " } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "iscx1:=(x,stepsiz)->x+step siz;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 107 "vecot:=array(1..25 1);ITER(de1,iscx1,iscy1,npoints=250,step=.01,x1=0,y1=0,datain=veco,nop lot,dataout=vecot);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 572 "Finally we can try to describe (mathematically!) the way the distance of the Brownian particle from the origin change s in time. To do this you need to superimpose a guess for the displace ment versus time against the output you have just obtained and to adju st your guess until it looks reasonable. The next line of commands jus t takes your output and makes a list of values of the displacement in \+ a form that can be plotted (as the plot structure p1). This is a piece of Maple for you to ignore! (But you need to edit the number of point s to be plotted and to enter it!) " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 105 "L:=[]:for i from 1 to 251 do l:=[i,vecot[i]]:L:=[[op(l)],op(L)] o d:p1:=plot(L,style=point):display(\{p1\});" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 363 "Exercise [6] The correct result f or the average displacement of a particle undergoing Brownian motion i s proportional to t^(1/2). Plot this, along with the ma gnitude of the displacements you obtained in p1, by editing the comman ds below. State your conclusion; in particuler, if necessary, explain why your result might not agree with the theory.\n" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "f(x):=x^.5:p2:=plot(f(x),x=0..250):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "display(\{p1,p2\});" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 119 "*******END OF UNIT 5**************************************\n******************* ****************************************" }}}}{MARK "0 1 0" 52 } {VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }