DECLARE SUB Pict2 () DECLARE SUB Pict3 () DECLARE SUB Pict4 () DECLARE SUB Pict5 () DECLARE SUB Pict1 () DECLARE SUB Pager () 86 PRINT CHR$(12); 88 REM 100 REM for further information about this program contact 200 REM conduit p.o. box 388 iowa city, iowa 52240 (319)353-5789 300 REM registry number: phy003 400 REM last revision 5/77 by neil s. ferguson for conduit 500 REM this program is distributed with the support of the national 600 REM science foundation grant no. Sed75-06596. Any opinion, findings, 700 REM conclusions, or recommendations expressed or implied are those 800 REM of the authors and do not necessarily reflect the views of the 900 REM national science foundation 1000 REM copyright by houghton-mifflin, 1976. All rights reserved. 1100 REM "Laplac" -- (basic program begins at line 290) 1200 REM 1300 REM 1400 REM description: 1500 REM 1600 REM this program solves laplace's equation over a 2d grid 1700 REM of points using simplest central difference quotients. 1800 REM 1900 REM 2000 REM 2100 REM instructions: 2200 REM 2300 REM this program is taken from "using computers in 2400 REM physics" by john r. merrill and should be used in conjunction 2500 REM with that text. The programs are distributed by project 2600 REM conduit. 2700 REM CALL Pict1 CALL Pager CALL Pict2 CALL Pager CALL Pict3 CALL Pager CALL Pict4 CALL Pager CALL Pict5 CALL Pager 12400 PRINT 12500 PRINT " v4=? ------------------------ V3=?" 12600 PRINT " 1 1" 12700 PRINT " 1 1" 12800 PRINT " 1 1" 12900 PRINT " 1 1" 13000 PRINT " 1 1 " 13100 PRINT " 1 1" 13200 PRINT " 1 1" 13300 PRINT " 1 1" 13400 PRINT " 1 1" 13500 PRINT " v1=? ------------------------- v2=?" 13600 CALL Pager 13700 PRINT 13800 PRINT "A you observed the diagram you should have noted that the" 13900 PRINT "potentials at v1(1,1), v2(10,1), v3(10,10) and v4(1,10) are" 14000 PRINT "not specified. This is to give you the option of choosing" 14100 PRINT "those values; for once those values are given they are" 14200 PRINT "connected by linear functions, such that the potentials" 14300 PRINT "on all four sides are defined. And since we are dealing" 14400 PRINT "with a 10 by 10 array we obtain an interior grid for cal-" 14500 PRINT "culating the interior potentials (for a deeper understanding" 14600 PRINT "of how the program works have a listing made at the end of" 14700 PRINT "this run)." 14800 PRINT 14900 PRINT "This is enough explanatory information, let's get on to" 15000 PRINT "the program." 15200 GOTO 15800 15700 STOP 15800 REM program begins 15900 PRINT 16000 PRINT "Do you wish to choose potentials v1,v2,v3 and v4 (yes or no)"; 16100 INPUT p$ 16200 IF p$ = "yes" THEN 16700 16300 IF p$ = "no" THEN 18300 16400 REM shouldn't be here 16500 PRINT "i don't understand your selection." 16600 GOTO 15900 16700 PRINT 16800 PRINT "What value do you choose for v1,v2,v3 and v4 (0 30 THEN 17900 17100 IF a < 0 THEN 17900 17200 IF b > 30 THEN 17900 17300 IF b < 0 THEN 17900 17400 IF c > 30 THEN 17900 17500 IF c < 0 THEN 17900 17600 IF d > 30 THEN 17900 17700 IF d < 0 THEN 17900 17800 GOTO 18900 17900 PRINT 18000 PRINT "This program is designed to be relatively rapid and simple; " 18100 PRINT "thus your selections must be between 0 and 30." 18200 GOTO 16700 18300 PRINT 18400 PRINT "I will v1=0, v2=10, v3=5 and v4=22." 18500 a = 0 18600 b = 10 18700 c = 5 18800 d = 22 CALL Pager 18900 DIM v(10, 10), u(10, 10) conv = .0005 19100 n = 0 19200 v1 = 0 19300 REM size of array 19400 i1 = 10 19500 j1 = 10 REM zero out v(i,j) FOR i = 1 TO i1 FOR j = 1 TO j1 v(i, j) = 0! NEXT j NEXT i 19600 REM v(i,j) boundry conditions 19700 v(1, 1) = a 19800 v(10, 1) = b 19900 v(10, 10) = c 20000 v(1, 10) = d 20100 REM slopes of straight lines connecting the corners 20200 m1 = (b - a) / (i1 - 1) 20300 m2 = (c - b) / (j1 - 1) 20400 m3 = (c - d) / (i1 - 1) 20500 m4 = (d - a) / (j1 - 1) 20600 REM calculate the potentials on the four sides 20700 REM side ab and dc 20800 FOR i = 2 TO i1 - 1 20900 v(i, 1) = v(1, 1) + m1 * (i - 1) 21000 v(i, 10) = v(1, 10) + m3 * (i - 1) 21100 v1 = v1 + v(i, 1) + v(i, 10) 21200 NEXT i 21300 REM side ad and bc 21400 FOR j = 2 TO j1 - 1 21500 v(1, j) = v(1, 1) + m4 * (j - 1) 21600 v(10, j) = v(10, 1) + m2 * (j - 1) 21700 v1 = v1 + v(1, j) + v(10, j) 21800 NEXT j 21900 v1 = v1 / (2 * (i1 + j1)) 22000 FOR i = 2 TO i1 - 1 22100 FOR j = 2 TO j1 - 1 22200 v(i, j) = v1 22300 NEXT j 22400 NEXT i 22500 v2 = 0 22600 FOR i = 1 TO i1 22700 FOR j = 1 TO j1 22800 u(i, j) = v(i, j) 22900 NEXT j 23000 NEXT i 23100 FOR i = 2 TO i1 - 1 23200 FOR j = 2 TO j1 - 1 23300 v(i, j) = (u(i + 1, j) + u(i - 1, j) + u(i, j + 1) + u(i, j - 1)) / 4 23400 IF v2 > ABS((v(i, j) - u(i, j)) / (v(i, j))) THEN 23600 23500 v2 = ABS((v(i, j) - u(i, j)) / v(i, j)) 23600 NEXT j 23700 NEXT i 23800 n = n + 1 23900 IF v2 > conv THEN 22500 24000 PRINT 24100 PRINT "# of iterations over grid ="; n PRINT "Convergence Criterion Fraction change is less than"; conv 24200 PRINT 24300 PRINT "The data is read in the following manner:" 24400 PRINT " v(1,10),v(2,10),. . .v(5,10),. . . . . . . v(10,10)" 24700 PRINT " v(1,9),. . . . . . . .v(5,9),v(6,9)=,. . . ." 25000 PRINT " v(1,8)=,. . . . ." 25100 PRINT " ." 25200 PRINT " ." 25300 PRINT " ." 25400 PRINT 25500 FOR j = j1 TO 1 STEP -1 25600 FOR i = 1 TO i1 25700 PRINT USING " ##.####"; v(i, j); 25800 NEXT i 25900 PRINT 26000 NEXT j PRINT " The PLOT.PLT file is being created, return to windows and run WPLOT." PRINT "You will see three 3D surface plots." REM now write out the plotting file OPEN "PLOT.PLT" FOR OUTPUT AS #2 REM example PRINT #2, "CLRDATA" PRINT #2, "CLRPLOT"; "" PRINT #2, "CLOSE" PRINT #2, PRINT #2, "TITLE 3-D Plot sliced "; j1; "x"; i1; " grid" PRINT #2, "ZLABEL Y" PRINT #2, "XLABEL X" PRINT #2, "YLABEL Potential in Volts " PRINT #2, "YSTART "; Vmin PRINT #2, "YEND "; Vmax PRINT #2, "YSTEP "; INT((Vmax - Vmin) / 20!) PRINT #2, "XSTART 1" PRINT #2, "XEND "; i1 PRINT #2, "XSTEP 1" PRINT #2, "XTICS 3" PRINT #2, "ZSTART 1" PRINT #2, "ZEND "; j1 PRINT #2, "ZSTEP 1" PRINT #2, "ZTICS 0" PRINT #2, "READZ "; i1 FOR k1 = 1 TO i1 PRINT #2, k1; NEXT k1 PRINT #2, PRINT #2, "READTAB "; i1; " "; i1 FOR j = 1 TO j1 STEP 1 PRINT #2, j; FOR i = 1 TO i1 PRINT #2, USING "##.#### "; v(i, j); NEXT i PRINT #2, NEXT j PRINT #2, "PLOT3D" PRINT #2, PRINT #2, "TITLE Rotated 3-Dimensional Plot" PRINT #2, "ROTATE" PRINT #2, "PLOT3D" PRINT #2, PRINT #2, "TITLE Surface Plot from LAPLACE" PRINT #2, "PLOT3D SURFACE" CLOSE #2 CALL Pager 26100 END SUB Pager 15302 PRINT , "Push Return to continue"; 15400 INPUT dum$ 15500 PRINT CHR$(12); END SUB SUB Pict1 2800 PRINT "Title: Laplace--a program for finding a solution to Laplace's" PRINT "equation." PRINT PRINT "The Laplacian operator, named after the mathematician Laplace, is" PRINT "defined as the" PRINT "divergence of the gradient of a scalar function (the student is" PRINT "assumed to have some background in vector analysis and electro-" PRINT "magnectic theory). It is a vector calculus second derivative and" PRINT "is a scalar. The function of which the laplacian is a second" PRINT "derivative is also a scalar, but the first derivative in be-" PRINT "tween is the gradient, a vector. By definition," PRINT PRINT " The Laplacian of f(x,y,z) is the div [grad(f)]" PRINT " which = Del^2 f(x,y,z)" PRINT "that is in Cartesian coordinates," PRINT PRINT " 2 2 2" PRINT " { d + d + d } fx,y,z) = 0" PRINT " ----- ------ ----- " PRINT " dx^2 dy^2 d z^2 " PRINT "This derivative is useful not only in the study of electro-" PRINT "magnetism, but also in the study of the propagation of all" PRINT "kinds of waves." END SUB SUB Pict2 4800 PRINT "Relative to the general electrostatic problem (which can" 4900 PRINT "be stated in terms of the laplacian), if you wish to find the" 5000 PRINT "electrostatic potential V throughout a charge-free region in" 5100 PRINT "space, you solve Laplace's equation with boundry conditions" 5200 PRINT "at the edges of the regions." 5300 PRINT 5400 PRINT "Question: Laplace's equation (for charge free regions) =" 5500 PRINT " del^2 V(x,y,z) = ? (If you are unable to input the" 5600 PRINT " correct answer, I can't permit you to continue)" 5700 INPUT p1 5800 IF p1 <> 0 THEN 6000 5900 GOTO 6500 6000 PRINT 6100 PRINT "This is unbelievable! I thought everyone would know the answer" 6200 PRINT "to that trivial question. Well, I'm generous--try it again," 6300 PRINT "as many times as you want?" 6400 GOTO 5300 6500 PRINT 6600 PRINT "I see you made it. Well let's continue." 6700 PRINT END SUB SUB Pict3 6800 PRINT "The general electrostatic problem can be stated in terms of" 6900 PRINT "poisson's equation (see program POISSON.BAS), Del^2 V= charge density/" 7000 PRINT "constant. This equation is simply a combination of Gauss's" 7100 PRINT "law and the fact that vector E=-del V(r). For charge-free" 7200 PRINT "regions, the equation reduces to laplace's equation," 7300 PRINT 7400 PRINT , "del^2 V(x,y,z)=0" 7500 PRINT 7600 PRINT "Now in most simple situations, you know the potential every-" 7700 PRINT "where around the surface of the region you are interested in;" 7800 PRINT "thus you use Poisson's or Laplace's equation to find V every-" 7900 PRINT "where. But, as you know, in general, the solution to partial" 8000 PRINT "differential equations is complicated whether you find the" 8100 PRINT "solution numerically or analytically. However, it happens that" 8200 PRINT "Laplace's equation and Poisson's equation, for many" END SUB SUB Pict4 8400 PRINT "situations is relatively simple and an adequate illustration" 8500 PRINT "of one way to treat partial differential equations on the" 8600 PRINT "computer." 8700 PRINT 8800 PRINT "The approach is: you approximate partial derivatives by " 8900 PRINT "difference quotients. Following are two simple (there are" 9000 PRINT "more complex ones) central difference approximations." 9100 PRINT 9200 PRINT "df/dx =(approx.) (f(x+c,y,z)-f(x-c,y,z))/2cx" 9300 PRINT "d**2f/dx**2 =(approx.) (f(x+c,y,z)-2f(x,y,z)+f(x-c,y,z))/(cx)**2" 9400 PRINT "(where c=delta increment and d=partial derivative)" 9500 PRINT 9600 PRINT "Using the difference equation you can find an expression for" 9700 PRINT "the value of one point in terms of values nearby. For sim-" 9800 PRINT "plicity, consider a two-dimensional region. Break that region" 9900 PRINT "up into a grid of points with equal grid spacing delta x and" 10000 PRINT "delta y. The value of V at a point (x,y) can be given in" END SUB SUB Pict5 10200 PRINT "terms of nearby grid-point values. An example: using the" 10300 PRINT "simple central difference approximation, Laplace's equation" 10400 PRINT "becomes (with delta x=delta y=d)," 10500 PRINT 10600 PRINT " (V(x+d,y)+v(x,y+d)+v(x-d,y)+v(x,y-d)-4v(x,y))/d**2=" 10700 PRINT " 0 (approx.)" 10800 PRINT 10900 PRINT " (V(x+d,y)+V(x,y+d)+V(x-d,y)+V(x,y-d))/4 = V(x,y)--approx." 11000 PRINT 11100 PRINT "This last equation says that the value of V at a point on" 11200 PRINT "the grid is approximately the average of the values at the four" 11300 PRINT "nearest neighbors." 11400 PRINT 11500 PRINT "The program Laplace illustrates one way to solve Laplace's" 11600 PRINT "equation in this manner. The potential on the edge of the" 11700 PRINT "rectangular region varies linearly. Initially, all the values" 11900 PRINT "of v inside the region are set to the average value on the" 12000 PRINT "boundry. The difference equation for v at each grid point is" 12100 PRINT "used over again at all interior points of the grid until the" 12200 PRINT "difference between iterations is small enough (say 0.1%). the bound" 12300 PRINT "points are always at the initially defined values." END SUB