Решение уравнений эллиптического типа

 

Program Direchle;

uses crt;

const n=900;

gr=3;

var z1,z2,poln,i_x,i_y,shag,x,y,ch,i,j,s:integer;

fi_le:text;

mas:array [-gr..gr,-gr..gr] of integer;

mas_height:array [-gr..gr,-gr..gr] of real;

sum:real;

begin

clrscr;

for j:=-gr to gr do {Задание граничных условий}

begin

writeln;

for i:=-gr to gr do

begin

if i=-gr then mas_height[j,i]:=0;

if j=gr then mas_height[j,i]:=0;

if j=-gr then

begin

mas_height[j,i]:=poln;

inc(poln);

end;

if i=gr then

begin

dec(poln);

mas_height[j,i]:=poln;

end;

end;

end;

writeln;

randomize;

assign(fi_le,'c:/direchle.dat');

rewrite(fi_le);

for i_x:=-gr+1 to gr-1 do

for i_y:=-gr+1 to gr-1 do

begin

ch:=0;

for i:=-gr to gr do

for j:=-gr to gr do

mas[i,j]:=0;

for i:=1 to n do {Расчет по блужданиям}

begin

x:=-i_x;y:=i_y;shag:=0;

while (x<>gr) and (y<>gr) and (x<>-gr) and (y<>-gr) do

begin

s:=random(9);

if (s=0) or (s=4) then begin inc(y);end;

if (s=1) or (s=5) then begin inc(x);end;

if (s=2) or (s=6) then begin dec(y);end;

if (s=3) or (s=7) then begin dec(x);end;

inc(ch);

inc(shag);

end;

mas[x,y]:=mas[x,y]+shag; {Количество попаданий}

end;

for z1:=-gr+1 to gr-1 do

sum:=sum+(mas[z1,gr]*mas_height[-z1,gr]/ch);

for z2:=gr-1 downto -gr+1 do

sum:=sum+(mas[gr,z2]*mas_height[-gr,z2]/ch);

mas_height[i_x,i_y]:=sum; {Значения искомой функции}

sum:=0;

end;

for j:=-gr to gr do {Формирование текстового файла}

begin

writeln(fi_le);

for i:=-gr to gr do

write(fi_le, mas_height[j,i]:2:3,' ');

end;

close(fi_le);

end.

Решение уравнений параболического типа

На примере уравнения теплопроводности

 

program montik;

uses crt;

const

stx=0.125;{шаг по длине}

stt=0.01; {шаг по времени}

mt=0.07; {время }

bts=9; {начальная температура на левом конце}

tle=0; { начальная температура на левом конце }

tpe=0; { начальная температура на левом конце }

Ni=500; {Число испытаний}

ld=0.5; {лямда}

var

mv:array[0..5] of real;

Mass:array[0..10,0..10] of real;

g:array[0..10,0..10] of integer;

N,M,Xx,Tt,x,t,i,vn,t1,x1:integer;

tv,rn,l:real;

f:text;

label Metka;

 

function Convert_St_Int(vh:real):integer;

var

vh_str:string;

err_vih,pos_s,st_int:integer;

begin

str(vh:2:2,vh_str);

pos_s:=pos('.',vh_str);

vh_str:=copy(vh_str,1,pos_s-1);

val(vh_str,st_int,err_vih);

Convert_st_int:=st_int;

end;

 

begin

clrscr;

assign(f,'c:\el.dat');

rewrite(f);

mv[0]:=5;

mv[1]:=ld*stt/(sqr(stx)+2*ld*stt);

mv[2]:=(stt-ld*stt)/(sqr(stx)+2*ld*stt);

mv[3]:=(sqr(stx)-2*stt+2*ld*stt)/(sqr(stx)+2*ld*stt);

mv[4]:=(stt-stt*ld)/(sqr(stx)+2*ld*stt);

mv[5]:=stt*ld/(sqr(stx)+2*ld*stt);

N:=Convert_St_Int(1/stx);

M:=Convert_St_Int(mt/stt);

Xx:=1;

Tt:=1;

randomize;

for x:=0 to N do

begin

Mass[x,0]:=bts;

end;

for t:=0 to M do

begin

Mass[0,t]:=tle;

Mass[N,t]:=tpe;

end;

for t:=1 to M do

begin {Цикл по времени}

for x:=1 to N-1 do

begin {Цикл по длине}

for i:=1 to Ni do

begin {цикл по испытаниям}

tv:=0;

rn:=random;

for vn:=1 to 5 do

begin {расчет по блужданиям}

if (tv<rn)and(rn<(tv+mv[vn])) then goto Metka;

tv:=tv+mv[vn];

end;{konec cikla po vn}

metka:

case vn of

1:Xx:=Xx-1;

2:

begin

Xx:=Xx-1;

Tt:=Tt-1;

end;

3:Tt:=Tt-1;

4:

begin

Xx:=Xx+1;

Tt:=Tt-1;

end;

5:Xx:=Xx+1;

end;{conec uslovii case}

if (Xx=0) or (Xx=N) or (Tt=0) then

begin

g[Xx,Tt]:=g[Xx,Tt]+1;

Xx:=x;

Tt:=t;

end;

end;{konec cicla po блужданиям}

for t1:=0 to t do

begin

Mass[x,t]:=Mass[x,t]+Mass[0,t1]*g[0,t1]/Ni;

Mass[x,t]:=Mass[x,t]+Mass[N,t1]*g[N,t1]/Ni

end;

for x1:=1 to N-1 do

begin

Mass[x,t]:=Mass[x,t]+Mass[x1,0]*g[x1,0]/Ni;

end;

for t1:=0 to 10 do

begin

for x1:=0 to 10 do

begin

g[x1,t1]:=0;

end;

end;

end;{conec cicla po X}

end;{conec cicla po T}

writeln;

l:=0;

write('X= |');

for x1:=0 to N do

begin

write(' ',l:2:2,' |');

l:=l+stx;

end;

writeln;

for x1:=0 to N do

begin

write('--------');

end;

writeln;

l:=stt*M;

for t:=M downto 0 do

begin

write('T=',l:2:2,' |');

l:=l-stt;

for x:=0 to N do

begin

write(' ',Mass[x,t]:2:2,' |');

write(f,Mass[x,t]:2:2,' ');

end;

writeln;

writeln(f);

end;

readln;

close(f);

end.

Метод Монте-Карло при моделировании

Задач нейтронной физики

Программа моделирования движения нейтронов через пластинку без учета сорта ядра.

 

Program berillii;

uses crt;

const a=7;

n=10000;

var mu,psi,x,e,bc,bt,p,h,g1,g2,g3,g4,z,z0:real;

i,j,ch_scvoz,ch_otr,ch_pogl,ch_rass:integer;

mu_disable,e_disable:boolean;

begin

clrscr;

randomize;

for j:=1 to 5 do

begin

ch_scvoz:=0;

ch_otr:=0;

ch_pogl:=0;

ch_rass:=0;

h:=3; g1:=1; z0:=0;

for i:=1 to n do

begin

g2:=random; g3:=random;

if not e_disable then e:=g1;

if not mu_disable then mu:=2*g3-1;

z:=z0-(ln(g2)*mu/e);

if z>h then

begin

inc(ch_scvoz); mu_disable:=false; e_disable:=false;

end;

if z<0 then

begin

inc(ch_otr); mu_disable:=false; e_disable:=false;

end;

if (z<=h) and (z>=0)

then

begin

if e<0.01 then bc:=4 else bc:=4-2.2*g1;

bt:=bc+(0.113/sqrt(g1));

p:=bc/bt;

g4:=random;

if g4<p then

begin

inc(ch_pogl);

mu_disable:=false;

e_disable:=false;

end

else

begin

x:=2*pi*g2;

psi:=(1+a*mu)/sqrt(1+2*a*mu+a*a);

e:=e*(a*a+2*a*mu+1)/(a+1)*(a+1);

mu:=mu*psi+cos(x)*sqrt((1-mu)*(1-psi*psi));

inc(ch_rass);

mu_disable:=true;

e_disable:=true;

end;

end;

g1:=random;

end;

writeln(ch_scvoz/n*100:5:0,');

writeln(ch_otr/n*100:5:0,'доля отраженных=');

writeln(ch_pogl/n*100:5:0, 'Доля поглощенных=’);

writeln(ch_rass/n*100:5:0,' доля проходящих');

end;

end.

Программа 2 моделирования прохождения нейтрона через вещество