Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Методичка по МС

.pdf
Скачиваний:
7
Добавлен:
22.03.2016
Размер:
865.98 Кб
Скачать

 

 

Компьютерное моделирование физических процессов

MatrixR^[1,3]:=0;

MatrixR^[2,3]:=0;

MatrixR^[3,3]:=1;

end; {SetMatrixR} {********************************************************************}

procedure ShowFigure(xx,yy:Word;p:Extended; c:Byte); {Вывод фигуры заданным цветом} {xx,yy - положение центра фигуры, p - угол поворота, c - цвет}

begin

SetColor(c); {Устанавливаем цвет для рисования} SetMatrixR(p); {Заполняем матрицу поворота на угол p} MatrixM^[1]:=Figure[1].X; {Заполняем матрицу первой точки} MatrixM^[2]:=Figure[1].Y;

MatrixM^[3]:=1;

MulMatrix(MatrixR,MatrixM,MatrixN); {Умножаем матрицу поворота на матрицу точки и получаем результат в матрице MatrixN} MoveTo(xx+Round(dl*MatrixN^[1]),yy+Round(dl*MatrixN^[2])); {Перемещаем перо в точку} for i:=2 to FN do {Для всех точек кроме первой}

begin

MatrixM^[1]:=Figure[i].X; {Заполняем матрицу точки} MatrixM^[2]:=Figure[i].Y;

MatrixM^[3]:=1;

MulMatrix(MatrixR,MatrixM,MatrixN); {Вычисляем положение точки после поворота} LineTo(xx+Round(dl*MatrixN^[1]),yy+Round(dl*MatrixN^[2]));{Линию в новую точку} end;

end; {ShowFigure} {********************************************************************}

begin

{Запускаем генератор случайных чисел}

Randomize;

ClrScr;

{Очистка экрана}

{Выделяем память под матрицы}

New(MatrixC);

New(MatrixA);

New(MatrixB);

New(MatrixM);

New(MatrixN);

New(MatrixR);

{----- Демонстрация работы с

матрицами -----}

{Замолняем MatrixA и MatrixB случайными числами}

for i:=1 to N do

 

for j:=1 to N do

 

begin

 

 

MatrixA^[i,j]:=10*Random;

 

MatrixB^[i,j]:=10*Random;

 

end;

{Выводим на экран матрицу MatrixA}

ShowSqrMatrix(MatrixA,2,3);

WriteLn;

{Выводим на экран матрицу MatrixB}

ShowSqrMatrix(MatrixB,2,3);

AddMatrix(MatrixA,MatrixB,MatrixC); {Складываем матрицы}

WriteLn;

{Выводим на экран результат суммирования}

ShowSqrMatrix(MatrixC,3,3);

ReadKey;

{Ожидаем нажатия любой клавиши}

{

----- Представление движения при помощи матриц -----}

gd:=Detect; gm:=0;

 

{Устанавливаем переменные, определяющие графический режим}

InitGraph(gd,gm,'');

{Устанавливаем графический режим}

x:=GetMaxX div 2; y:=GetMaxY div 2;

{Координаты центра фигуры}

repeat

 

{Повторяем до нажатия любой клавиши}

 

ShowFigure(x,y,phi,0);

{Выводим фигуру черным цветом (стираем)}

 

phi:=phi+0.1;

 

{Изменяем угол}

 

ShowFigure(x,y,phi,15);

{Выводим фигуру белым цветом}

 

Delay(500);

 

{Задержка}

 

until keypressed;

 

{}

 

ReadKey;

 

{Ожидаем нажатия любой клавиши}

CloseGraph;

 

{Закрываем графику}

{Освобождаем память, выделенную под матрицы}

Dispose(MatrixA);

Dispose(MatrixB);

Dispose(MatrixC);

Dispose(MatrixM);

Dispose(MatrixN);

Dispose(MatrixR);

end.

 

 

 

41

Компьютерное моделирование физических процессов

Приложение 3. Движение тела брошенного под углом к горизонту

Program

Movement;

{Подключаемые модули}

uses

CRT,Graph;

const

g=9.81;

 

{Ускорение свободного падения}

var

dt=1E-4;

: Integer;

{Квант времени}

gd,gm

{Переменные для инициализации графики}

 

x,y

: Extended;

{Координаты тела}

 

vx,vy,v : Extended;

{Проекции и модель скорости тела}

begin

alpha

: Extended;

{Угол к горизонту, под которым брошено тело}

 

ReadLn(v);

{Вводим начальную скорость}

Write('v=');

Write('alpha='); ReadLn(alpha);{вводим угол}

x:=0; y:=0;

 

{Начальные координаты}

vx:=v*cos(alpha);

{Вычисляем проекции на оси координат}

vy:=v*sin(alpha);

 

gd:=Detect; gm:=0;

{Определяем тип графического режима}

InitGraph(gd,gm,'');

{Устанавливаем графический режим}

repeat

 

 

{Вычисляем горизонтальную скорость через dt}

vx:=vx+0;

 

vy:=vy-g*dt;

 

{Вычисляем вертикальную скорость через dt}

x:=x+vx*dt;

 

{Координаты теле через dt}

y:=y+vy*dt;

 

 

PutPixel(Round(x),GetMaxY-Round(y),15);{Ставим точку на экране}

until (keypressed)or(y<0);

{Выход из цикла по нажатию на клавишу или при падении

тела на поверхность Земли}

 

ReadKey;

 

 

{Ожидание нажатия любой клавиши}

CloseGraph;

 

{Закрытие графики и выход из программы}

end.

 

 

 

42

Компьютерное моделирование физических процессов

Приложение 4. Колебания математического маятника

Program

Oscillation;

{Подключаемые модули}

uses

CRT,Graph,DOS;

const

R=200;

 

{Длина нити = радиус обращения}

var

gd,gm

: Integer; {Переменные для инициализации графики}

 

Save1C

: Pointer; {Переменная для хранения адреса прерывания 1C}

 

Phi

: Extended;{Угол отклонения}

 

PhiDot

: Extended;{Угловая скорость}

 

Phi2Dot : Extended;{Угловое ускорение}

 

Omega

: Extended;{Частота колебаний маятника}

 

Delta

: Extended;{Параметр затухания}

{$F+}

t,dt

: Extended;{Время и квант времени}

 

 

{Команда компилятору для организации дальних вызовов}

procedure MyTimer; interrupt;{Процедура, вызываемая по прерыванию от таймера}

var dPhiDot : Extended;

{Изменение угловой скорости}

procedure ShowM(c:Byte);

{Процедура вывода маятника на экран}

var

xx,yy : Integer;

{}

 

begin

 

 

{Цвет вывода маятника}

SetColor(c);

 

MoveTo(GetMaxX div 2,GetMaxY div 2);{Переходим в средину экрана}

xx:=Round(R*cos(3*pi/2+phi));

{Вычисляем смещение по оси x}

yy:=Round(R*sin(3*pi/2+phi));

{Смещение по оси y}

LineRel(xx,-yy);

 

{Рисуем линию}

Circle(GetMaxX div 2+xx,GetMaxY div 2-yy,4);{Вывод окружности - тело}

end;

 

 

 

 

begin

 

 

{Отобразить маятник черным - стереть}

ShowM(0);

 

dPhiDot:=-(2*Delta*PhiDot+Sqr(Omega)*sin(Phi))*dt;{Изменение угловой скорости}

Phi2Dot:=dPhiDot/dt;

 

{Угловое ускорение}

Phi:=Phi+(PhiDot+dPhiDot/2)*dt;

{Угол поворота}

PhiDot:=PhiDot+dPhiDot;

 

{Угловая скорость}

t:=t+dt;

 

 

{Время}

ShowM(15);

 

 

{Вывод белым цветом}

end;

 

 

 

{Команда компилятору для организации ближних вызо-

{$F-}

 

 

 

вов}

 

 

 

 

begin

 

ReadLn(Phi);

{Начальный угол}

Write('Phi0=');

Write('PhiDot='); ReadLn(PhiDot);{Начальная скорость}

Write('Omega=');

ReadLn(Omega); {Частота колебаний}

Write('Delta=');

ReadLn(Delta); {Параметр затухания}

dt:=1/18.2;

 

 

{Квант времени, определяется частотой таймера}

gd:=Detect; gm:=0;

 

{Параметры графического режима}

InitGraph(gd,gm,'');

 

{Устанавливаем графику}

GetIntVec($1C,Save1C);

 

{Получаем вектор прерывания 1С и сохраняем его в

переменной Sava1C}

 

{Устанавливаем вектор прерывания 1С на процедуру

SetIntVec($1C,@MyTimer);

 

MyTimer}

 

 

 

 

Repeat Until KeyPressed;

 

{Пустой цикл до нажатия на любую клавишу}

SetIntVec($1C,Save1C);

 

{Восстанавливаем вектор прерывания}

CloseGraph;

 

 

{Закрываем графику и выходим из программы}

end.

 

 

 

 

43

Компьютерное моделирование физических процессов

Приложение 5. Линии напряженности электростатического поля

program

LineE;

 

 

{Подключаемые подули}

uses

Crt,Graph;

 

type

TPoints=Record

{Запись для хранения данных о положении и величине заряда}

 

X,Y,q : Extended;

const

end;

 

{Количество зарядов}

NZ=3;

 

 

k=9E9;

 

{Электрическая постоянная}

 

dl=1E-1;

{Длина вектора dl}

 

TE=1E1;

{Расстояние, в пределах которого строятся силовые линии}

 

R1=3;

 

{Радиус, с которого строятся силовые линии}

 

RMin=0.1; {Расстояние от заряда, на котором оканчивается линия напряженности}

 

dPhi=2*pi/20; {Шаг изменения угла, с которым строятся линии}

var

gd,gm : Integer;

{Переменные инициализации графики}

 

PosZ

: Array[1..NZ] of TPoints; {Массив зарядов}

 

ii

: Integer;

{Переменная цикла}

MaxX,MaxY : Integer;

{Размер экрана}

{**********************************************************************}

procedure

ShowAllLines(w:Integer); {Построение силовых линий от заряда w}

var Ex,Ey,E

: Extended;

{Проекции и модуль вектора напряженности}

x,y

: Extended;

{Координаты точки линии напряженности}

i

: Integer;

{Переменная цикла}

r

: Extended;

{Расстояние между точками}

dx,dy

: Extended;

{Проекции вектора dl на оси координат}

phi

: Extended;

{Угол, определяющий начало силовой линии}

b,b1

: Boolean;

{Переменные условий}

key

: Char;

{Переменная для хранения кода нажатой клавиши}

begin

{Начальный угол, под которым строится первая линия}

Phi:=0.0001;

b:=True;

{Переменная определеяющая переход к следующей линии}

repeat

 

 

{----- Определяем координаты начала новой линии -----}

if b then

{Если построение линии закончено}

begin

x:=PosZ[w].X+R1*cos(Phi); {Координаты начала новой линии}

y:=PosZ[w].Y+R1*sin(phi);

{Линия не закончена}

b:=False;

b1:=False;

{Угол для построения следующей линии}

phi:=phi+dphi;

end;

 

{----- Строим линию напряженности -----}

Ex:=0; Ey:=0;

{Обнуляем переменные}

for i:=1 to NZ do

{Перебираем все заряды}

begin

{Выполняем действия для i-ого заряда}

r:=sqrt(sqr(x-PosZ[i].X)+sqr(y-PosZ[i].Y)); {Расстояние от заряда до точки}

E:=k*PosZ[i].q/sqr(r);

{Модуль напряженности в точке}

Ex:=Ex+E*(x-PosZ[i].X)/r;

{Проекции на оси координат}

Ey:=Ey+E*(y-PosZ[i].Y)/r;

 

if r<RMin then b1:=True; {Если расстояние до заряда меньше Rmin, то линия закончилась}

end;

dx:=Ex/sqrt(sqr(Ex)+sqr(Ey))*dl;

{Приращение координаты x}

dy:=Ey/sqrt(sqr(Ex)+sqr(Ey))*dl;

{Приращение координаты y}

MoveTo(Round(x),Round(y));

{Строим линию в новую точку}

LineRel(Round(dx),Round(dy));

{Координаты новой точки}

x:=x+dx; y:=y+dy;

if (b1)or(x>MaxX+TE)or(x<-TE)or(y<-TE)or(y>MaxX+TE)then b:=True;

{Если линия окончилась на заряде или ушла за пределы рассматриваемой области, то переход к следующей линии}

until (keypressed)or(Phi>2*pi+dPhi);

44

{Рамка} {Запуск генератора случайных чисел}

Компьютерное моделирование физических процессов

{Выход если нажата любая клавиша или линия окончилась на заряде или вышла за пределы рассматриваемой области}

end;

begin

gd:=Detect; gm:=0; InitGraph(gd,gm,''); {Инициализация графики} MaxX:=GetMaxX; {Размер экрана в пикселях} MaxY:=GetMaxY;

Rectangle(0,0,GetMaxX,GetMaxY);

Randomize;

for ii:=1 to NZ do {Перебираем все заряды} begin

PosZ[ii].X:=Random*GetMaxX; {Случайным образом выбираем координаты}

PosZ[ii].Y:=Random*GetMaxY;

{И величину заряда}

PosZ[ii].q:=Random-0.7;

if PosZ[ii].q>0 then SetColor(15) else SetColor(14);

Circle(Round(PosZ[ii].X),Round(PosZ[ii].Y),3); {Выводим на экран заряды разным цветом}

end;

 

for ii:=1 to NZ do

{Перебираем все заряды}

begin

{Определяем цвет}

SetColor(8+ii);

if PosZ[ii].q>0 then ShowAllLines(ii); {Строим силовые линии от всех положитель-

ных зарядов}

end;

 

ReadKey;

{Ожидаем нажатия любой клавиши}

CloseGraph;

{Закрываем графический режим и выходим из программы}

end.

 

45

Компьютерное моделирование физических процессов

Приложение 6. Распространение света в среде с градиентом АПП

program

GradienOptic;

{Подключаемые модули}

uses

CRT,Graph;

const

MaxN=3;

{Максимальное значение АПП}

 

d=1E-2;

{Расстояние между лучами, см. рис. 7.5}

 

dt=1E-11;

{Квант времени}

 

c=3E8;

{Скорость света в вакууме}

var gd,gm

: Integer {Переменные для инициализации графики}

i,j

: Integer;

{Переменные циклов}

phi

: Extended;

{Угол распространения луча}

dphi

: Extended;

{Изменение угла распространения луча}

R1,R2,x1,y1,x2,y2 : Extended; {Радиусы волн и координаты точек}

{********************************************************************}

{Функция вычисления COS угла с использованием математического }

{сопроцессора, позволяет значительно ускорить вычисления } function FastCos(c:Extended): Extended; assembler;

asm

fld

tbyte ptr [bp+4]

db 0D9h,

0FFh

end;

{********************************************************************} {Функция вычисления SIN угла с использованием мат. сопроцессора}

function

FastSin(c:Extended): Extended; assembler;

asm

tbyte ptr [bp+4]

fld

db 0D9h, 0FEh end;

{********************************************************************}

function fn(_x,_y:Extended):Extended;{Возвращает значение АПП в точке} begin

fn:=2+SQR(cos(_y/100+pi/2)); {Различные варианты изменения АПП}

{fn:=2+cos(_y/100+pi/2);}

{fn:=1+_y/768*2}

{if sqr(500-_x)+sqr(350-_y)>100000 then fn:=1 else fn:=3;}

{if (_x<100) or(_x>300) then fn:=1 else fn:=1.6;}

end;

{********************************************************************}

procedure ShowLines; {Вывод световых лучей} begin

x2:=x1+d*FastSin(phi); {Вычисление координат второго луча} y2:=y1-d*FastCos(phi);

repeat

{Вычисление радиусов сферических волн}

R1:=c/fn(x1,y1)*dt;

R2:=c/fn(x2,y2)*dt;

 

dPhi:=(R1-R2)/d;

{Изменение угла распространения света}

Phi:=Phi-dPhi;

{Новый угол распространения света}

x1:=x1+R1*FastCos(Phi); {Новые координаты точек первого} y1:=y1+R1*FastSin(Phi);

x2:=x2+R2*FastCos(Phi); {и второго лучей} y2:=y2+R2*FastSin(Phi);

PutPixel(Round(x1),Round(y1),100); {Рисование только первого луча}

46

Компьютерное моделирование физических процессов

until (keypressed)or(x1>GetMaxX)or(y1>GetMaxY)or(x1<0)or(y1<0); {Вычисления завершаются если нажата любая клавиша или луч вышел} {за пределы экрана}

end;

{********************************************************************}

begin

gd:=InstallUserDriver('Vesa256',nil); {Драйвер VESA256.BGI}

gm:=3;

{Установка

{Разрешение 1024х768}

InitGraph(gd,gm,'');

графического режима}

{Используется VESA-драйвер, позволяющий

использовать 256 цветов}

{и высокое разрешение экрана}

 

 

for i:=0 to 63 do SetRGBPalette(i,i,i,i);

{Цвета с номурами от 0 до 63 устанавливаются оттенками серого} {Белый цвет 63, черный 0}

SetRGBPalette(100,0,0,255); {Цвет 100 - Синий}

SetRGBPalette(110,255,0,0); {Цвет 110 - Красный}

for i:=1 to GetMaxX+1 do

{Перебираем все точки экрана от левого}

for j:=1 to GetMaxY+1 do

{верхнего угла до правого нижнего}

PutPixel(i,j,Round(63-(fn(i,j)-1)*63/MaxN)); {Ставим точки, цвет которых определяет АПП в данной точке. Среда с большим АПП - темнее}

for i:=1 to 20 do

{Строим 20 лучей}

begin

 

 

x1:=1; y1:=200+i*20; {С координатами x1, y1 идущих под углом phi}

phi:=30*pi/180;

 

ShowLines;

 

end;

 

 

ReadKey;

{Ожидание нажатия любой клавиши}

CloseGraph; {Закрытие графического режима и выход из программы} end.

47

Компьютерное моделирование физических процессов

Приложение 7. Явление радиоактивного распада

Program

RadioActiv;

{Подключаемые модули}

uses

CRT,Graph;

const

N0=10000;

{Начальное количество нераспавшихся ядер}

 

Lambda=0.01;

{Постоянная распада}

var

dt=1E-2;

{Квант времени}

gd,gm

: Integer;

{Переменные для установки графики}

 

N

: LongInt;

{Количество не распавшихся ядер}

 

N1

: LongInt;

{Количество распавшихся ядер}

 

i

: LongInt;

{Переменная цикла}

begin

t

: Extended;

{Время}

 

 

{Запускаем генератор случайных чисел}

Randomize;

 

gd:=Detect; gm:=0;

{Определяем тип графического режима}

InitGraph(gd,gm,'');

{Устанавливаем графический режим}

N:=N0;

 

 

{Начальное количество нераспавшихся ядер}

repeat

 

 

{Число распадов равно нулю}

N1:=0;

 

 

for i:=1 to N do

{В цикле по количеству нераспавшихся ядер}

if Lambda*dt>Random then inc(N1); {Проверяем, произошел ли распад, если да, то увеличиваем счетчик распавшихся ядер}

{Случайность распада определяется при сравнении со случайным числом}

N:=N-N1;

{Число нераспавшихся ядер уменьшается на число распадов}

t:=t+dt;

{Увеличиваем время}

PutPixel(Round(t),Round(GetMaxY*(1-N/N0)),15); {Ставим точку на графике}

until (KeyPressed)or(N<=0);

{Выход по нажатию на клавишу или по распаду всех}

ReadKey;

{Ожидание нажатия любой клавиши}

CloseGraph;

{Закрыть графику и выйти}

end.

 

48

Компьютерное моделирование физических процессов

Приложение 8. Темы и план лабораторных занятий

Тема 1. Матричное представление движения.

1.Создание модуля для выполнения основных операций с матрицами.

1.1.Создайте новый тип для квадратной матрицы размером 3 на 3.

1.2.Создайте новый тип для матрицы–столбца размером 1 на 3.

1.3.Процедуры и функции для выполнения основных операций над матрицами — сложение и умножение квадратных матриц, умножение квадратной матрицы на матрицу–столбец.

1.4.Создание матриц основных операций — матрицы параллельного переноса, поворота, деформации, зеркального отражения.

2.Проверка работоспособности модуля и правильности проводимых вычислений

2.1.Напишите программу, демонстрирующую сложное движение плоской фигуры. Движение организуется при помощи операций с матрицами.

2.2.Расширить возможности модуля, добавив функции работы с матрицами 4 на 4

и1 на 4. Разработать матрицы операций для пространственного движения тела.

Тема 2. Интегрирование дифференциальных уравнений.

1.Движение тела, брошенного под углом к горизонту.

1.1.Напишите программу, моделирующую движение тела, брошенного под углом к горизонту без учета сил вязкого трения. Проверьте правильность работы по уравнениям известным из кинематики.

1.2.Измените программу, учтя в ней силу вязкого трения.

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

отображены время подъема, время падения, время движения, максимальная

высота подъема, дальность полета, конечная скорость (угол между вектором скорости и горизонтальной осью, модуль скорости), скорость тела в верней

точке траектории. Все величины должны быть определены численно, а не

формулам кинематики.

2.Колебания математического маятника.

2.1.Напишите программу, моделирующую колебания математического маятника с учетом сил вязкого трения и внешней вынуждающей силы.

2.2.Постройте уравнения зависимости угла поворота, угловой скорости и углового ускорения от времени.

2.3.Пронаблюдайте явления резонанса и биений.

2.4.Измените программу так, чтобы колебания маятника происходили в реальном

времени (используя прерывания от системного таймера). Тема 3. Моделирование движения тел Солнечной системы.

1.Напишите программу, моделирующую движение тел Солнечной системы.

2.Пронаблюдаете движение тел Солнечной системы.

3.Используя модель докажите справедливость законов Кеплера, определите сино-

дический и сидерический период обращения планет, значения больших и малых полуосей орбит, секторные скорости в перигелии и афелии.

4.Используя модель, определите амплитуду колебаний Солнца.

Тема 4. Построение эквипотенциальных и силовых линий электростатического поля.

49

Компьютерное моделирование физических процессов

1.Напишите программу, демонстрирующую эквипотенциальные линии электростатического поля.

1.1.Используя генератор случайных чисел, расставьте на экране N зарядов, величину каждого заряда также задайте случайным числом, заряды должны быть обоих знаков.

1.2.Рассчитайте потенциал каждой точки поля (за исключением областей в непосредственной близости от зарядов).

1.3.Постройте карту потенциального поля, используя оттенки серого. 1.4.Постройте эквипотенциальные линии в оттенках красного (синего) цвета.

2.В разработанную ранее программу добавьте функцию построения силовых линий 2.1.Область вблизи каждого положительного заряда разбейте на сектора и через каждый сектор постройте силовую линию. При разбиении помните условие

из электростатики, что густота силовых линий определяет величину напря-

женности поля.

2.2.Проводите построение линий до тех пор, пока они не закончатся на отрица-

тельных зарядах.

2.3.Пронаблюдаете перпендикулярность линий напряженности и силовых линий.

2.4.Постройте стандартные картины силовых полей: два одинаковых одноименных

или разноименных заряда. Сравните вид полей с представленным в физике. Тема 5. Моделирование распространения света в среде с градиентом АПП.

1.Задайте функцией закон изменение АПП в пространстве, учтите, что АПП должен

быть больше 1. Наглядно отобразите среду на экране в оттенках серого цвета.

2.Напишите процедуру, моделирующую распространение света в заданной среде.

3.Пронаблюдайте ход светового луча в различных средах. 3.1.Среда с плавным изменением АПП.

3.2.Среда с резкой границей АПП (закон преломления света). 3.3.Среды с несколькими резкими границами АПП: линзы, призмы.

Тема 6. Создание модели ядерного реактора.

1.Напишите программу, реализующую ядерный реактор.

2.Пообдерите параметры для ручного управления ядерным реактором.

3.Модернизируйте программу для автоматического управления реактором и полу-

чения различных мощностей.

4.Определите КПД полученного реактора.

5.Смоделируйте аварийную ситуацию с нарушением функций управления. Оцени-

те энергетический выход неуправляемой цепной реакции.

50