ПРИЛОЖЕНИЕ 3. РАСПЕЧАТКИ ПРОГРАММ НА ЯЗЫКЕ ПАСКАЛЬ: ЧАСТЬ 2

Ниже приводятся версии перевода восьми программ из ч. 2 данной книги с языка True BASIC иа Паскаль. Программы пропускались с использованием компилятора Macintosh Pascal. Для других версий Паскаля необходимо внести минимальные изменения. Применяемые в программах на Паскале графические процедуры входят в состав стандартного программного обеспечения компьютера Macintosh и могут вызываться из других языков программирования, реализованных на этом компьютере. Краткое описание графических процедур для компьютера Macintosh приведено в табл. Д1.

ГЛАВА 10

program integ (input, output);

(* вычисляется интеграл от f(x) в пределах от х = а до х = b *) var

a,  b,  h, area : real; n : integer;

procedure initial (var a, b, h : real;

var n : integer);

begin

a := 0.0;           (* нижний предел *)

b := 0.5 * 3.14159;      (* верхний предел *)

write('число интервалов = '); readln(n);

h :=~(b - а) / п (* шаг интегрирования *)

end;

function f (x : real) : real; begin

f := cos(x) end;

procedure rectangle (a, b, h : real;

n : integer;

var area : real);

var

x, sum : real;

i : integer; begin x := a; sum := 0.0; for i := 0 to n - 1 do begin

sum := sum + f(x); x := x + h end;

area := sum * h end;

procedure outdat (area : real); begin

writeln('площадь = ', area : 12 : 7) end;

begin   (* основная программа *)

initial(a, b,  h, n);

rectangle(a, b, h, n, area);

outdat( area) end.

ГЛАВА 11

program random_walk (input, output);

{* моделирование Монте-Карло одномерного случайного блуждания *) const

xmax = 64; type

list = array[-64. .64] of real;

var

prob : list;

p : real;            (* вероятность шага вправо *)

x, xcum, x2cum, itrial, ntrial, nstep : integer;

procedure initial (var р : real;

var nstep,  ntrial, xcum, x2cum : integer; var prob : list);

var

i : integer; begin

ntrial := 100;    (* число испытаний *)

p := 0.5;           {* вероятность шага вправо *)

write(' максимальное число шагов = ');

readln( nstep);

write('начальное случайное число (отрицательное целое) = '); read ln( iseed);

(* начальное обнуление данных *) xcum := 0; x2c0m := 0;

for i := -nstep to nstep do prob[i] := 0

end;

function rnd : real;

(* выдается случайное число от 0 до 1 с использованием средств *) (* компьютера Macintosh. Функция random выдает случайное число *) (* от -32768 до 32767. Более хорошие генераторы случайных чисел *) (* можно найти в книге Пресса У. и др.. Численные рецепты *) begin

rnd := (random + 32768.0)/{32768.0 + 32767.0) end;

procedure walk (var x : integer;

p : real;

nstep : integer);

var

istep : integer; begin

x := 0;  (* начальное положение для каждого испытания *)

for istep := 1 to nstep do if rnd <= p then x := x + 1

else

x := x - 1

end;

procedure data (var x, xcum, x2cum : integer; var prob : list);

begin

xcum := xcum + x;

x2cum := x2cum + x * x;

prob[x] := prob[x] + 1 end;

procedure average (nstep, ntrial, xcum, x2cum : integer; var prob : list);

(• средние значения для nstep-шагового блуждания *) var

norm, xbar, x2bar, sigma, variance, probability : real; x : integer; begin

norm := 1.0 / ntrial; xbar := xcum • norm; x2bar := x2cum • norm; (or x := -nstep to nstep do begin

probability := prob[x] • norm; writeln(x : 6, probability : 10 : 5) end;

variance := x2bar - xbar * xbar; sigma := sqrt( variance);

miteln('среднее смещение = ', xbar : 10 : 5); writeln('sigma   = ', sigma : 10 : 5) end;

begin   (* основная программа *)

initial(p, nstep, ntrial, xcum, x2cum, prob); for itrial := 1 to ntrial do begin

walk(x, p, nstep);

data(x, xcum, x2cum, prob) («сбор данных после nstep шагов*) end;

average( nstep, ntrial, xcum, x2cum, prob) end.

ГЛАВА 12

program site_perc (input, output);

(* рисуются конфигурации ячеечной перколяции *)

const

boxlength = 180;

marginy = 10;

marginx = 100; type

list = array[1..50, 1..50] of real;

var

L : integer; r : list;

procedure initial (var L : integer); begin

write('линейный размер решетки = '); readln(L);

(* рисование ячеек решетки *)

framerect( margin, margin, boxlength + margin, boxlength + margin) end;

function rnd : real;

{* выдается случайное число от 0 до 1  с помощью функции *) (* random, которая выдает случайное число от -32768 до 32767 *) (• Более хорошие генераторы случайных чисел можно найти в *) (• книге Пресса У. и др.. Численные рецепты *) begin

rnd := (random + 32768.0)/(32768.0 + 32767.0) end;

procedure lattice (L : integer;

var r :  I ist);

(* каждой ячейке придается случайное число и рисуется решетка *) var

col,  row, х, у : integer; scale : real; begin

scale := boxlength / L;

for row := 1 to L do    (* рисование ячеек решетки *)begin

у := round((row - 0.5) * scale) + marginy; (* с каждой ячейкой связывается квадрат со стороной 1*) (or col := 1 (о L do begin

x := round(scale * (col - 0.5)) + marginx;

(* каждой ячейке придается случайное число *)

r[col, row] := rnd;

{* в каждой ячейке рисуется точка *) moveto(x, у); lineto(x, у)

end

end

end;

procedure configuration (L : integer;

var г :  I ist); (• занятие ячеек с данной вероятностью р *) var

s : list;

р, scale : real;

row, col, x, у, size : integer; begin

for row := 1 to L do . for col := 1 to L do s[row, col] := 0; scale := boxlength / L; while p <= 1.0 do begin

write(' вероятность p = '); readln( p);

(* половина стороны квадрата для занятых ячеек *) size := round{0.5 • scale); for row := 1 to L do begin

у := round(scale * (row - 0.5)) + marginy; for col := 1 to L do

if (r[col, row] < p) and (s[col, row] <> 1.0) then

begin   (• новые занятые ячейки *)

х := round(scale * (col - 0.5)) + marginx;paintrect(y - size, x - size, у + size' x + size);s[col, row] := 1.0      (• ячейку заняли *)

end

end

end

end;

begin   (* основная программа *)

initial(L);

(* каждой ячейке придается случайное число и рисуется решетка *) lattice(L, г);

configuration L, г)      (* занятие ячеек с данной р *)

end.

ГЛАВА 13

program single_cluster (input, output);

(* вычисление фрактальной размерности перколяционного кластера *)

(* с помощью алгоритма Хаммерсли-Лиса-Александровича *)

type

list = array[1..200] of integer;

var

L : integer; p : real; s : list;

procedure parameter (var L : integer;

var p : reaI);

begin

write('максимальный радиус кластера = '); readln( L);

write('вероятность занятия ячейки = '); readln( p) end;

function rnd : real;

(• вычисление случайного числа от 0 до 1 *) begin

rnd := (random + 32768.0) * (1.5259е-5) end;

procedure grow (L : integer;

p : rea I;

var s :   I ist); (* расчет перколяционного кластера *) var

i> k. x, y, xn, yn, nper, nn, r : integer; r2 : real;

px, py : array[1..5000] of integer;

lat : array[-50..50, -50..50] of integer;

nx, ny : arrayf_1..4] of integer; (* в массивах nx и ny указаны положения ближайших соседей *) (• относительно каждой ячейки *) begin

for г := 1 to 2*L do s[r] := 0;

for x := -L to L do for у := -L to L do lat[x,y] := 0; lat[0, 0] := 1; nx[1] := 1; nx[2] := -1; nx[3] := 0; nx[4] := 0; пуП] := 0; ny[2] := 0; пуГЗ] := 1; ny[4] := -1; for i := 1 to 4 do begin

px[i] := nx[i]; py[i] := пуП] end;

nper := 4;         (* начальное число ячеек периметра •)

repeat

к := 1 + trunc(rnd * nper); («случайный выбор ячейки периметра*) х := рх[к]; у := руГ>]; if (rnd < р) then begin

(* site tested and occupied *) lat[x, y] := 1; r2:=x*x + y*y;

r := round(sqrt(r2));    (* расстояние от начала координат *) s[r] := s[r] + 1;        {* число ячеек иа расстоянии г *) (• вновь занятая ячейка заменяется в массивах рх и ру *) (• последней ячейкой в списке ячеек периметра *) рх[к] := рх[прег]; PYfk] := prfnper]; nper := nper - 1;

(• отыскание новых ячеек периметра *) for nn := 1 to 4 do begin

xn := x + nx[nn]; yn := у + nv[nn]; if ((latfxn, yn] = 0) and (abs(xn) <= L) and (abs(yn) <= L)) then begin

nper := nper + 1; px[nper] := xn; PYfnper] := yn end

end

end else begin

(* ячейка проверяется, но не занимается *) lat[x, у] := -1; px[k] := pxfnper]; PyM := PYlnper]; nper := nper - t end

until (nper < 1)   (* пока не будут проверен весь периметр *)

id;

procedure plotg (L : integer;

s : list);

const    (* минимальные используемые экранные координаты *)

yO = 190; xO = 50;

var

xmax, ymax, sea lex, scaley, hashx, hashy : integer; x, y, r, xcoor, ycoor, mass : integer; begin

mass := s[1];

xmax := trunc(ln(L) + 1);

ymax := trunc(ln(4 * L * L) + 1);

scalex := 300 div xmax;

scaley := 190 div ymax;

hashy -.= 10;   (* размеры делений *)

hashx := 15;

(* вычерчивание осей *) moveto(x0, yO - ymax * scaley); Iineto(x0, yO);

Iineto(x0 + xmax * scalex, yO); (* маркировка осей *) moveto(250, 270); drawstring{'ln{r)'); moveto(3, 80); drawstring('ln(M)'); (* нанесение делений на осях *) for х := 1 to xmax do begin

xcoor := x * scalex + xO; moveto( xcoor, y0); lineto( xcoor, yO - hashy) end;

for у := 1 to ymax do begin

ycoor := yO - scaley * y; moveto{ xO, ycoor); Iineto(x0 + hashx, ycoor) end;

{* рисование графика 1п{заключенной массы) от 1п(г) *) for г := 2 to L do begin

mass := mass + s[r];

x := xO + trunc(scalex * ln(r));

Y := trunc(yO - sea ley * In(mass));

paintrect(Y - 1, x - 1, у + 1, x + 1)

end

end;

begin   (* основная программа *)

parameter( L, p);

grow{L, p, s);

plot(L, s) end.

ГЛАВА 14

program entropy (input, output);

{* вычисляется энтропия методом подсчета совпадений Ma *) type

list = аггауТ_1..10] of integer; state = array[0. .2000] of integer;

var

micro : state;

right,  left : list;

nl, nr, nexch : integer;

procedure start (var n I, nr, nexch : integer;

var left,  right : list;

var micro : state); (* ввод параметров и выбор начальной конфигурации частиц *) var

ir, il, N : integer; begin

write{' полное число частиц = '); readln(N);

writej'число частиц слева = '); readln(nl);

nr := N - nl;     (* число частиц с правой стороны *)

micro[0] := 0; for il := 1 to nl do begin

left[il] : = il;     (* список номеров частиц с левой стороны *) micro[0] := micro[0] * 2 + 2   («начальное микросостояние*) end;

for ir := 1 to nr do

right[ir] : = ir + nl;    {* список частиц с правой стороны *) write('число обменов = '); readln( nexch) end;

function rnd : real;

(* формируется случайное число в интервале от 0 до 1 *) begin

rnd := (random + 32768.0) / (32768.0 + 32767.0) end;

function power (n -. integer) : integer;

{* выдается 2 в степени n *)

var

i, p : integer; begin

P := 1;

for i := I to n do

p := p * 2; power := p end;

procedure exch (var n I, nr, nexch : integer;

var left,  right : list;

var micro : state); (• обмен частицы слева с номером ileft *) («с частицей справа с номером iright *) var

iexch, ileft, iright, jleft, jright : integer; begin

for iexch := 1 to nexch do begin

(• индексы массивов выбираем случайным образом *) ileft := trunc(rnd * nl + 1); iright := trunc(rnd * nr + 1); jleft := left[ileft]; jright := right[iright];

left[ileft] := jright;   (* новый номер частицы слева *) right[iright] := jleft; {* новый номер частицы справа *) (• определяем новое микросостояние •) micro[iexch] := micro[iexch - 1] + power{jright); micro[iexch] :- micro[iexch] - power(jleft) end

end;

procedure outdat (nexch : integer;

var micro : state); (* вычисляются частота совпадений и энтропия *) (• полное число сравнений •) var

iexch, jexch, ncomp, ncoin : integer; rate, S : real; begin

ncomp := nexch • (nexch - 1) div 2; ncoin := 0;

(* сравниваем микросостояния •) for iexch := 1 to nexch - 1 do

for jexch := iexch + 1 to nexch do

if {microfiexch] = microfjexch]) thenncoin := ncoin + 1;rate := ncoin / ncomp;        (» частота совпадений »)

if (rate > 0) then

S := ln{1.0 / rate); writeln('оценка энтропии = ', S) end;

begin   (* основная программа •)

start(nl, nr, nexch, left, right, micro);

exch(nl, nr, nexcl, left, right, micro);  {» обмен частиц *)

outdat( nexch, micro)     (• расчет частоты совпадений и энтропии *) end.

ГЛАВА 15

program idealdemon (input, output);

(* алгоритм демона для одномерного идеального классического газа *) type

list = array[1..100] of real;

var

vel : list;

N, nmcs, imcs, i : integer;

esystem, edemon, vtot, dvmax, escum, vcum, edcum, accept : real;

procedure initial (var N, nmcs : integer;

var esystem, edemon, vtot, dvmax : real; var vel  : list);

var

i : integer; vinitial : real; begin

write{' число частиц = '); readln(N);

write{'число шагов Монте-Карло - '); readln( nmcs);

write(' начальная энергия системы = '); readln( esystem);

edemon := 0;   {* энергия демона *)

write(' максимальное изменение скорости = '); readln( dvmax);

vinitial := sqr{2 * esystem / N);

{» у всех частиц одинаковые начальные скорости •) vtot := 0.0; for i := 1 to N do begin

vel[i] := vinitial;

vtot := vtot + vinitial   {* полная скорость системы *)

end

end;

function rnd : real; begin

rnd := (random + 32768.0) / (32768.0 + 32767.0) end;

procedure changes (N : integer;

var esystem, edemon, vtot, dvmax, accept : real; var vel : list);

var

dv, vtrial, de : real; ip : integer;

begin

dv := (2 * rnd - 1) * dvmax;       {* пробное изменение скорости *)ip := frunc(rnd » N + 1);      {* берем случайную частицу »)

vtrial := vel[ip] + dv;   (* пробная скорость *)

de := 0.5 * (vtrial * vtrial - vel[ip] * vel[ip]); («изменение энергии*) if (de <= edemon) or (de <= 0) then begin

velfj'p] := vtrial;

vtot := vtot + dv;        (* полная скорость системы •)

accept := accept + 1; edemon := edemon - de; esystem := esystem + de end

end;

procedure data (var esystem, edemon, vtot, escum, vcum, edcum : real); begin

edcum      edcum + edemon; escum := escum + esystem; vcum := vcum + vtot end;

procedure averages (N, nmcs integer;

escum, vcum, edcum, accept : real);

var

norm, edave, vave, esave : real; begin

norm := 1.0 / (nmcs » N);

edave := edcum • norm;          (• средняя энергия демона •)

accept := accept * norm;         {* коэффициент принятия *)

(• средние системы на частицу *) norm := norm / N;

esave := escum * norm;    («средняя энергия на частицу системы*) vave := vcum » norm;     («средняя скорость на частицу системы *) writeln{' средняя энергия демона = ', edave : 10 : 5); writeln('средняя энергия системы на частицу = ', esave : 10 : 5); writeln('средняя скорость на частицу = ', vave : 10 : 5); writeln('коэффициент принятия = ', accept : 10 : 5) end;

begin   {* основная программа *)

initial(N, nmcs, esystem, edemon, vtot, dvmax, vel); accept := 0.0; escum := 0.0; vcum := 0.0; edcum := 0.0; for imcs := 1 to nmcs do for i := 1 to N do begin

changes(Nr esystem, edemon, vtot, dvmax, accept, vel); (* накапливаем данные после каждого испытания *) data( esystem,  edemon, vtot,  escum, vcum, edcum) end;

averages(N, nmcs, escum, vcum, edcum, accept) end.

ГЛАВА 16

program Ising2 (input, output);

(• алгоритм Метрополиса для двумерной модели Изинга *) type

lattice = array[1..32, 1..32] of integer; list = array[-4. .4] of real;

var

spin : lattice; w : list;

N, L, imcs, nmcs : integer;

T, E, M, ratio, ecum, e2cum, mcum, m2cum : real; function rnd : real;

(* формируется случайное число в интервале от 0 до 1 •) begin

rnd := (random + 32768.0) / (32768.0 + 32767.0) end;

procedure initial (var N, L, nmcs : integer;

var spin : I attice; var T, E, M : real; var w : list);

i, j, up, right, sum : integer; e4, e8 : real; begin

write('линейный размер решетки = '); readln( L);

N := L * L;      (* число спинов «)

ууп1е('число шагов Монте-Карло на спин = '); readln( nmcs);

write{'приведенная температура = '); readln(T);

(* случайная начальная конфигурация •) for j := 1 to L do for i := 1 to L do begin

if (rnd < 0.5) then spin[i, j] := 1

else

spin[i,  j] := -1; M      M + spin[i, j]     (* суммарная намагниченность *) end;

(• вычисление начальной энергии Е *) for j := 1 to L do begin

if (j = L) then up := 1

else

up := j + 1; for i := 1 to L do begin

if (i = L) then right := 1

else

right := i + 1; sum := spinfi,  up] + spin[right, j]; E := E - spin[i, j] * sum       (• полная энергия •)

end

end;

(• вычисление вероятностей перехода *)

(• индекс массива w равен сумме спинов четырех соседей *)

е4 := ехр{-4.0 / Т); е8 := е4 * е4; w[4] := е8; w[-4] := е8; w[2] := е4; w[-2] := е4 end;

procedure pbc (var i, j, L, sum : integer;

var spin : latt i ce); (* находится сумма значений соседних спинов; •) (• используются периодические краевые условия *) var

left, right, up, down : integer; begin

if (i = 1) then

left := spin[L, j]

else

left := spin[i - 1, j]; if (i = L) then

right := spin[1, j]

else

right := spin[i + 1, j]; if (j = 1) then

down := spin[i, L]

else

down := spin[i,  j - 1]; if (j = L) then

up := spin[i, 1]

else

up :- spin[i, j + 1]; sum := left + right + up + down

end;

procedure accept (i, j, sum : integer;

var spin :   I at t i ce; var ratio, M, E : real);

begin

spin[i, j] := -spin[i, j]; ratio := ratio + 1; M := M + 2 * spin[i, j]; E := E - 2 * spin[i, j] * sum end;

procedure Metropolis (N, L :  i nteger;

var spin : lattice; var E, M, ratio: real; var w : list);

var

ispin, i, j, sum : integer; begin

for ispin := 1 to N do begin

(* случайный разыгрыш координат спнна *) i := trunc(L * rnd + 1); j := trunc(L * rnd + 1);

(* находим значения ближайших соседних спинов •)

(* периодические краевые условия •)

pbc(i, j, L, sum, spin);

if (spin[i, j] * sum <= 0) then

accept(i, j, sum, spin, ratio, M, E) else if (rnd < w[sum]) then

accept(i, j, sum, spin, ratio, M, E)

end

end;

procedure data (var E, M, ecum, e2cum, mcum, m2cum : real);

{* накапливание данных после каждого шага Монте-Карло на спин •)

begin

ecum := ecum + Е;

e2cum := e2cum + E * E;

mcum := mcum + M;

m2cum := m2cum + M » M end;

procedure outdat (N, nmcs : integer;

ecum, e2cum, mcum, m2cum, ratio : real);

var

norm, eave, e2ave, save, s2ave : real; begin

norm := 1.0 / (nmcs » N); (* средние на спин *) ratio := ratio * norm; eave := ecum » norm; e2ave := e2cum * norm; save : = mcum • norm; s2ave := m2cum * norm; writeln('коэффициент принятия := ', ratio); writeln('средняя энергия на спин := ', eave); writeln(' средний квадрат энергии на спин := ', e2ave); writeln(" средняя намагниченность := ', save); writeln(' средний квадрат намагниченности := ', s2ave) end;

begin   (* основная программа *)

initial(N, L, nmcs, spin, T, E, M, w); (or imcs := 1 to nmcs do begin

Metropolis(N, L, spin, E, M, ratio, w); data(E, M, ecum, e2cum, mcum, m2cum) end;

outdat(N, nmcs, ecum, e2cum, mcum, m2cum, ratio) end.

ГЛАВА 17

program eigen (input, output); var

VO, a, xmax, dx, xO, yO, xscale, \tca\e : real;

procedure parameter (var VO, a, xmax, dx : real); begin

write(' глубина ямы = '); readln(VO);

write{' полуширина ямы = ');

readln( а);

write{'максимальное выводимое на график значение х = '); readln( xmax);

write(' размер шага dx = '); readln( dx); end;

procedure plot_potential (VO, a, xmax : real;

var xscale, yscale, *0» yO : real); (* задание масштаба графика н построение потенциала*) begin

xscale := 512 / (2.0 » xmax); yscale := 300 / (2.0 * 10.0); xO := xscale • xmax; y0 := yscale * Ю.0;

moveto(trunc(xO - xmax • xscale, trunc(y0 + 6 * yscale)); lineto{trunc(x0 - a * xscale, frunc(y0 + 6 * yscale)); lineto(trunc(x0 - a * xscale, frunc(y0 + 9 » yscale)); lineto(frunc{xO + a * xscale, trunc(y0 + 9 * yscale)); Iinefo(frunc(x0 + a * xscale, trunc(y0 + 6 » yscale)); Iineto(trunc(x0 + xmax • xscale, trunc(y0 + 6 * yscale)); end;

function V (x, V0, a : real) : real;

(• вычисляется функция потенциала V(x) *)

begin

if x > a then

V         := V0

else

V         :- 0.0

end;

procedure Euler (V0, a, dx, xmax, xscale, yscale, xO, yO : real); var

i, j, parity : integer; x, phi, dphi, d2phi, E .- real; begin

write(' четная или нечетная функция (1 или -1) '); readln( parity);

repeat

write('E = '); readln(E);

(• задание начальных значений при х = 0 *) if parity = -1 then begin

phi := O.O;

dphi := 0.0

end

else

begin

phi := 1.0;

dphi := 0.0 end; x := 0.0;

repeat  (* вычисление волновой функции •)

x := x + dx;

d2phi := 2.0 * phi * (V(x, VO, a) - E);

dphi := dphi + d2phi * dx;

phi := phi + dphi * dx;

(• plot wave function •)

i := trunc(x0 + xscale • x);

j := trunc(y0 + yscale * phi);

moveto(i, j);

lineto(i,j);

i := trunc(x0 - xscale • x);j := trunc(y0 + yscale * phi * parity);moveto{ i, j);lineto( i, j);until (x > xmax) or (abs(phi) > 10);until E = 0 (• для выхода из программы введите Е = 0 *)

end;

begin   (• основная программа •)

parameter( V0, a, xmax, dx);

plot_potential(V0, a, xmax, xscale, yscale, xO, yO); Euler(V0, a, dx, xmax, xscale, yscale, xO, yO) end.