Основы компьютерного моделирования

Размещено на http://www.allbest.ru/

ОСНОВЫ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ

Учебное пособие

Глазов: ГГПИ, 2005



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



1. Принципы компьютерного моделирования

2. Методы численного интегрирования и дифференцирования

3. Моделирование систем с одной степенью свободы

4. Модель двумерного движения материальной точки

5. Модель движения системы материальных точек

6. Моделирование колебаний связанных осцилляторов

7. Моделирование распространения одномерной волны

8. Моделирование явлений переноса

9. Моделирование автоволновых процессов



1. ПРИНЦИПЫ КОМПЬЮТЕРНОГО МОДЕЛИРОВАНИЯ

Компьютерное моделирование является одним из эффективных методов изучения физических систем. Часто компьютерные модели проще и удобнее исследовать, они позволяют проводить вычислительные эксперименты, реальная постановка которых затруднена или может дать непредсказуемый результат. Логичность и формализованность компьютерных моделей позволяет выявить основные факторы, определяющие свойства изучаемых объектов, исследовать отклик физической системы на изменения ее параметров и начальных условий.

Компьютерное моделирование требует абстрагирования от конкретной природы явлений, построения сначала качественной, а затем и количественной модели. За этим следует проведение серии вычислительных экспериментов на компьютере, интерпретация результатов, сопоставление результатов моделирования с поведением исследуемого объекта, последующее уточнение модели и т.д.

К основным этапам компьютерного моделирования относятся: постановка задачи, определение объекта моделирования; разработка концептуальной модели, выявление основных элементов системы и элементарных актов взаимодействия; формализация, то есть переход к математической модели; создание алгоритма и написание программы; планирование и проведение компьютерных экспериментов; анализ и интерпретация результатов.

Различают аналитическое и имитационное моделирование. Аналитическими называются модели реального объекта, использующие алгебраические, дифференциальные и другие уравнения, а также предусматривающие осуществление однозначной вычислительной процедуры, приводящей к их точному решению. Имитационными называются математические модели, воспроизводящие алгоритм функционирования исследуемой системы путем последовательного выполнения большого количества элементарных операций.

Принципы моделирования состоят в следующем [3]:

1. Принцип информационной достаточности. При полном отсутствии информации об объекте построить модель невозможно. При наличии полной информации моделирование лишено смысла. Существует уровень информационной достаточности, при достижении которого может быть построена модель системы.

2. Принцип осуществимости. Создаваемая модель должна обеспечивать достижение поставленной цели исследования за конечное время.

3. Принцип множественности моделей. Любая конкретная модель отражает лишь некоторые стороны реальной системы. Для полного исследования необходимо построить ряд моделей исследуемого процесса, причем каждая последующая модель должна уточнять предыдущую.

4. Принцип системности. Исследуемая система представима в виде совокупности взаимодействующих друг с другом подсистем, которые моделируются стандартными математическими методами. При этом свойства системы не являются суммой свойств ее элементов.

5. Принцип параметризации. Некоторые подсистемы моделируемой системы могут быть охарактеризованы единственным параметром: вектором, матрицей, графиком, формулой.

Компьютерное моделирование систем часто требует решения дифференциальных уравнений [1-10]. Важным методом является метод сеток, включающий в себя метод конечных разностей Эйлера. Он состоит в том, что область непрерывного изменения одного или нескольких аргументов заменяют конечным множеством узлов, образующих одномерную или многомерную сетку, и работают с функцией дискретного аргумента, что позволяет приближенно вычислить производные и интегралы. При этом бесконечно малые приращения функции и приращения ее аргументов заменяются малыми, но конечными разностями.

2. МЕТОДЫ ЧИСЛЕННОГО ИНТЕГРИРОВАНИЯ И ДИФФЕРЕНЦИРОВАНИЯ

1. Задача. Для известной функции , определите первую и вторую производные в точке с координатой , а также интеграл в интервале от до

2. Теория. Пусть задана функция . Разобьём интервал от до на элементарные отрезки длиной

,

получив конечное множество узлов сетки

где а - число узлов. В результате функция непрерывного аргумента будет заменена функцией дискретного аргумента . Тогда левая, правая и центральная разностные производные первого порядка в точке с координатой соответственно равны:

Чем меньше шаг сетки , тем выше точность найденных производных.

Вторая производная определяется из выражения:

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

,

длина каждой из которых равна

то площадь будет примерно равна:

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

Более точный метод трапеций заключается в том, что каждая -ая полоска рассматривается как трапеция высотой с длинами оснований

,

поэтому ее площадь равна

Интеграл функции равен сумме всех элементарных площадей этих трапециевидных полосок:

Метод Монте-Карло нахождения площади криволинейной трапеции под кривой состоит в следующем. Представим себе прямоугольник, ограниченный пределами интегрирования и , осью и горизонталью внутри которого находится эта криволинейная трапеция. Площадь прямоугольника равна . Задавая случайным образом координаты поместим внутрь прямоугольника точек. Подсчитаем число точек, оказавшихся внутри криволинейной трапеции, то есть удовлетворяющих условию

Площадь криволинейной трапеции будет во столько раз меньше площади выбранного прямоугольника, во сколько раз меньше . Поэтому при дробь стремится к пределу, равному искомому интегралу:

3. Компьютерная программа. Самостоятельно составьте алгоритмы нахождения производных и интегралов. Ниже представлены примеры программ. Первая программа позволяет вычислить первую и вторую производные функции

в точке с координатой , а также найти ее интеграл в интервале от до методом трапеций. Вторая программа определяет интеграл функции в интервале от 0 до 1 методом Монте-Карло.

program PROGRAMMA1_1;

uses crt;

var x,y1,y2,y3,a,b,h,S :real;

Function Funct(x:real):real;

begin {Задание функции}

Funct:=x*x*x-x*x+3;

end;

BEGIN {Основная программа}

clrscr; x:=3; h:=0.001;

y1:=Funct(x-h);

y2:=Funct(x);

y3:=Funct(x+h);

Writeln('Первая производная ', (y2-y1)/h:3:3);

Writeln('Вторая производная ', (y1-2*y2+y3)/(h*h):3:3);

a:=1; b:=3; x:=a; S:=0;

Repeat {Интеграл}

S:=S+0.5*(Funct(x)+Funct(x+h))*h; x:=x+h;

until x>b;

Writeln('Интеграл ',S:3:3);

Repeat until KeyPressed;

END.

program PROGRAMMA1_2;

uses crt;

const NN=10000;

var x,y,xx,yy: real;

n,i: integer;

function Funct(x:real):real;

begin Funct:=x*x; end;

BEGIN {Основная программа}

clrscr; Randomize; n:=0;

for i:=1 to NN do

begin

x:=Random(1000)/1000;

yy:=Random(1000)/1000;

if yy<Funct(x) then n:=n+1;

end;

writeln('Интеграл равен ',n/NN);

Repeat until KeyPressed;

END.

4. Задания для студентов.

1. Точка движется прямолинейно со скоростью

Численными методами и аналитически определите ускорение точки в момент времени 1 с и пройденный путь за время от 1 до 3 с. Повторите расчеты при различных значениях шага и сравните результаты.

2. Имеется пластинка толщиной ограниченная кривой

и прямой Ее плотность есть функция координаты

:

где - произвольный коэффициент пропорциональности. Определите ее площадь и массу методом Монте-Карло.

3. Определите координаты центра масс пластины толщины , ограниченной прямыми

,

плотность которой равна

4. Пластина толщиной имеет форму круга радиуса Ее плотность с ростом расстояния до центра убывает по закону

.

Методом численного интегрирования определите момент инерции пластины относительно оси проходящей через ее центр и лежащей в ее плоскости.

8. Постройте кривую зависимости излучательной способности абсолютно черного тела от частоты при постоянной температуре выражаемую формулой Планка:

где и - постоянные коэффициенты. Постройте график при различных .

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

10. Методом численного интегрирования найдите интегральную светимость абсолютно черного тела, взяв интеграл

Выполните расчеты для различных температур и убедитесь в справедливости закона Стефана-Больцмана: Интегральная светимость абсолютно черного тела пропорциональна четвертой степени его температуры:

3. МОДЕЛИРОВАНИЕ СИСТЕМ С ОДНОЙ СТЕПЕНЬЮ СВОБОДЫ

1. Задача. Имеется физическая система с одной степенью свободы, состоящая из инерционного элемента массой , упругого элемента жесткостью и диссипативного элемента с коэффициентом сопротивления . Определить отклик системы , а также ее первую и вторую производные на внешнее воздействие

если известны начальные условия

2. Теория. Из второго закона Ньютона следует линейное неоднородное дифференциальное уравнение второго порядка:

Процессы, происходящие в последовательном колебательном контуре, состоящем из последовательно соединенных резистора , конденсатора и катушки индуктивности , на который подано напряжение , описываются уравнением

где - заряд, проходящий по цепи,

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

,

происходят процессы, описываемые уравнением:

где - напряжение на резисторе, емкости и индуктивности, - производная силы тока по времени. При этом реализуется электромеханическая аналогия "сила-ток".

Механические величины	Аналогия сила-напряжение	Аналогия сила-ток
Сила	Напряжение	Ток
Смещение	Заряд	Напряжение
Скорость	Сила тока	Производная напряжения
Масса	Индуктивность	Емкость
Жесткость	Величина обратная емкости	Величина обратная индуктивности
Коэффициент сопротивления	Сопротивление	Сопротивление

Итак, неоднородное диффуравнение второго порядка описывает широкий класс задач, изучаемых в школьном и вузовском курсе физики: движение связанного с пружиной тела в вязкой среде под действием произвольной силы , протекание тока через последовательно или параллельно соединенные резистор, конденсатор и катушку индуктивности, подключенные к источнику ЭДС или источнику тока .

Характер движения механической системы зависит от действующей на нее внешней силы. При этом могут быть рассмотрены следующие ситуации:

· внешняя сила отсутствует;

· внешняя сила постоянна;

· внешняя сила изменяется по гармоническому закону;

· внешняя сила изменяется по произвольному периодическому закону;

· внешняя сила изменяется по произвольному непериодическому закону.

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

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

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

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

Рассмотрим алгоритм численного решения уравнения.

1. Задают параметры физической системы зависимость внешнего воздействия от времени , а также начальные условия и шаг по времени

2. Начало цикла по Дают приращение по времени: переменной присваивают значение

3. Определяют ускорение, скорость и координату тела в момент

4. Результаты вычислений выводят на экран в числовом виде либо строят соответствующие точки на координатной плоскости.

5. Возвращение к операции 2. Если цикл по закончился, - выход из цикла.

4. Компьютерная программа. При запуске программа рисует графики зависимостей координаты

проекции скорости

и проекции ускорения

.

Некоторые строчки программы заключены в скобки "(*" и "*)". Убрав скобки и активизировав соответствующие операторы, можно промоделировать различные явления.

program PROGRAMMA2;

uses dos, crt, graph;

Const Fm=10;w=5;m=2;r=0;k=0;

Mx=20; Mv=40; Ma=8; Mf=2; Mt=100;

dt=0.00006;

Var x,v,a,F,t : Real;

j,xx,vv,aa,FF,tt,Gd,Gm : Integer;

BEGIN

Gd:= Detect;

InitGraph(Gd, Gm, 'c:\bp\bgi');

if GraphResult<>grOk then Halt(1);

t:=0; v:=0; x:=-3;

line(30,300,650,300);

line(31,500,31,10);

OutTextXY(50,20,'X, V, A');

Repeat

begin {Задание функции F=F(t)}

t:=t+dt; (* F:=Fm*sin(w*t); *)

(*If sin(w*t)<0 then F:=0;

If sin(w*t)>0 then F:=Fm;*)

F:=0; If t<1 then F:= Fm;

If t>3 then F:=-Fm;

a:=(F-r*v-k*x)/m; x:=x+v*dt; v:=v+a*dt; tt:=round(t*Mt);

xx:=round(x*Mx); vv:=round(v*Mv); aa:=round(a*Ma); FF:=round(F*Mf);

circle(30+tt,300-xx,1); circle(30+tt,300-vv,1); circle(30+tt,300-aa,2);

end;

until KeyPressed;

CloseGraph;

END.

5. Задания для студентов.

1. На точку массы действует скачкообразно изменяющаяся сила

Если

Если

Если

Исследуйте движение точки, проанализируйте получившиеся графики зависимостей

2. Промоделируйте движение материальной точки, движущейся в вязкой среде под действием постоянной силы, направленной вдоль оси

:

при начальных условиях

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

3. Создайте модель переходного процесса в цепи, содержащей резистор и катушку индуктивности подключенные к источнику постоянного напряжения, при условии, что . Исследуйте аналогичный переходный процесс в цепи, содержащей последовательно соединенные резистор и конденсатор.

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

.

Убедитесь в том, что ускорение изменяется в противофазе с координатой, а скорость опережает координату на причем амплитуды колебаний уменьшаются по экспоненте. Проведите серию вычислительных экспериментов при различных начальных условиях системы.

5. Промоделируйте движение осциллятора в случае сильного затухания при

Убедитесь, что в этом случае движение будет апериодическим.

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

7. Промоделируйте работу сглаживающего -фильтра при подаче на него пульсирующего напряжения, получающегося в результате однополупериодного выпрямления. Для этого необходимо задать следующие параметры колебательной системы кг, Нс/м, Н/м, рад и закон изменения вынуждающей силы:

Если

Если

Где

- период колебаний. Убедитесь в том, что с ростом индуктивности уменьшается коэффициент пульсаций тока и напряжения на резисторе. Изучите зависимость амплитуды пульсаций от индуктивности , сопротивления нагрузки и частоты импульсов

8. Изучите работу интегрирующей цепи, состоящей из последовательно соединенных резистора и конденсатора, с которого снимается выходное напряжение. Видно, что при подаче на цепь прямоугольных импульсов

Если

Если

заряд конденсатора, а значит и напряжение на нем, возрастает пропорционально интегралу от входного напряжения.

Так как в программе осуществляется деление на (аналог индуктивности ), то значение этого параметра должно быть очень малым, но не равным нулю.

9. Промоделируйте движение тела в вязкой среде (), на которое в момент времени начинает действовать внешняя гармоническая сила

Эта ситуация соответствует переходному процессу, происходящему при подключении активно-индуктивной нагрузки к источнику переменного напряжения. При переходный ток стремится к принужденному току, изменяющемуся с той же частотой, что и приложенная ЭДС и отстающему от нее на некоторую фазу.

10. Создайте программу, моделирующую процессы, происходящие в колебательной системе в случае, если на нее действует периодически изменяющаяся сила, частота которой пропорциональна времени:

,

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

.

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



4. МОДЕЛЬ ДВУМЕРНОГО ДВИЖЕНИЯ МАТЕРИАЛЬНОЙ ТОЧКИ

1. Задача. Материальная точка массой движется в силовом поле

,

при этом на нее действует сила вязкого трения

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

моделирование точка одномерный автоволновый

2. Теория. Примерами подобного движения являются движение точки, в однородном силовом поле, в центральном силовом поле сил притяжения или отталкивания, в центральном поле сил упругости и т.д. При этом могут быть учтены силы вязкого трения.

Проанализируем основные ситуации.

1. Движение в однородном поле. Во всех точках пространства вектор силы имеет постоянные проекции на оси координат. При отсутствии силы трения точка движется по параболе, а при ее наличии - по более сложной кривой.

2. Движение в центрально-симметричном поле, действующем по закону обратных квадратов. На точку с координатами действует сила

Ее проекции на оси координат:

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

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

,

лежащая в плоскости экрана и направленная перпендикулярно вектору скорости. Введем угол , который образует вектор с осью . Проекции силы Лоренца на координатные оси:

Заряженная частица описывает окружность. При наличии тормозящей силы радиус окружности уменьшается.

4. Движение частицы в скрещенных электрическом и магнитном полях. Пусть силовые линии электрического поля лежат в плоскости экрана и направлены вверх, а силовые линии магнитного поля направлены к нам перпендикулярно вектору .

Если заряд частицы положительный, то на него со стороны электрического поля действует постоянная сила, направленная вверх. Чтобы учесть ее влияние необходимо к вертикальной проекции силы Лоренца прибавить постоянное слагаемое :

Если начальная скорость частицы равна нулю, то траекторией ее движения является циклоида.

3. Алгоритм. Пусть в момент времени материальная точка имеет координаты и проекции скорости Запишем второй закон Ньютона:

Отсюда следует, что проекции ускорения точки в момент времени равны:

Определив координаты и проекции скорости точки в момент времени , можно повторить процедуру вычисления требуемое количество раз и построить траекторию движения точки.

Построим алгоритм модели.

1. Задают массу материальной точки , коэффициент вязкости , начальные координаты и проекции скорости силовое поле

,

а также шаг по времени

2. Начало цикла по Дают приращение по времени: переменной присваивают значение

3. Определяют ускорение, скорость и координату тела в следующий момент времени:

4. Результаты вычислений выводят на экран в числовом виде либо строят соответствующие точки на координатной плоскости.

5. Возвращение к операции 2. Если цикл по закончился, -- выход из цикла.

4. Компьютерная программа. Предлагаемая компьютерная программа позволяет изучить движение материальной точки в различных силовых полях с учетом действующей на точку силы трения.

program PROGRAMMA3;

uses crt, graph;

var v, B, q, F, Fx, Fy : real;

r, x, y, vx,vy,ax,ay : real; Gd, Gm, i: integer;

const M=500; mm=100; dt=0.005; rr=0.1; k=2;

Begin

Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');

if GraphResult <> grOk then Halt(1);

line(320,240,640,240); line(320,240,320,0); circle(320,240,5);

x:=100; y:=120; vx:=1; vy:=-2;

Repeat

begin

{--Заданние силового поля--}

(* Fy:=3; Fx:=0; *)

(* Fx:=-k*x; Fy:=-k*y; *)

(* r:=sqrt(x*x+y*y); F:=M*mm/(r*r);

Fx:=-F*x/r; Fy:=-F*y/r; *)

B:=2; q:=1; F:=B*v*q; v:=sqrt(vx*vx+vy*vy);

Fx:=F*vy/v; Fy:=-F*vx/v;

(* B:=2; q:=1; F:=B*v*q; v:=sqrt(vx*vx+vy*vy);

Fx:=F*vy/v; Fy:=-0.5-F*vx/v; *)

{--Расчет скоростей и ускорений--}

ax:=(Fx-rr*vx)/mm; ay:=(Fy-rr*vy)/mm;

vx:=vx+ax*dt; vy:=vy+ay*dt; x:=x+vx*dt; y:=y+vy*dt;

circle(round(x)+320,240-round(y),2); setcolor(12);

circle(round(x)+320,240-round(y),1); setcolor(15);

end;

until KeyPressed;

CloseGraph;

END.

5. Задания для студентов.

1. Материальная точка массы брошена под углом к горизонту в однородном поле тяжести. Изучите траекторию ее движения при отсутствии силы вязкого трения и при ее наличии.

2. Убедитесь в том, что время подъема материальной точки, движущейся под действием силы тяжести в вязкой среде, меньше времени спуска. Воспользуйтесь тем, что в наивысшей точке подъема проекция скорости на ось меняет свой знак на противоположный.

3. Промоделируйте движение точки в поле центральных сил упругости

в случае, когда на точку действует сила вязкого трения и когда она равна нулю. По какой траектории движется точка?

4. Исследуйте движение точки в поле сил притяжения, действующих по закону обратных квадратов

Промоделируйте ситуации, в которых точка движется по гиперболе, параболе, эллипсу. Изучите характер движения искусственного спутника Земли, входящего в верхние слои атмосферы, на который действует сила вязкого трения.

5. Изучите движение точки в поле сил отталкивания. Промоделируйте опыт Резерфорда по отклонению альфа-частиц ядрами атомов золота. Меняя прицельный параметр, проведите серию компьютерных экспериментов.

6. Промоделируйте движение заряженной частицы в камере Вильсона, помещенной в однородное магнитное поле. Учтите, что по мере своего движения частица теряет кинетическую энергию. Магнитное поле направлено перпендикулярно плоскости экрана.

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

8. Изучите движение материальной точки в гравитационном поле двух массивных тел. Проведите компьютерные эксперименты при заданных начальных координатах и скорости точки.

9. Задайте произвольное силовое поле

и промоделируйте движение частицы в нем.

5. МОДЕЛЬ ДВИЖЕНИЯ СИСТЕМЫ МАТЕРИАЛЬНЫХ ТОЧЕК

1. Задача. Имеется система материальных точек с массами , взаимодействующих друг с другом с внутренними силами

на каждую из которых действует внешняя сила

.

Исходя из начальных радиус-векторов и скоростей , определите координаты и скорости материальных точек в последующие моменты времени.

2. Теория. Рассмотрим механическую систему из материальных точек. Основной закон динамики:

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

Дифференциальное уравнение второго порядка может быть представлено двумя дифференциальными уравнениями первого порядка. Имеем:

Зная внешние и внутренние силы, действующие на каждую материальную точку, можно определить их ускорения. Исходя из координат и скоростей точки в момент времени можно рассчитать координаты и скорости точки в следующий момент времени

Так как любая механическая система - совокупность взаимодействующих между собой материальных точек, то эта модель чрезвычайно широка и охватывает большое число механических систем. Помимо моделирования одномерного и двумерного движения материальной точки, данный метод позволяет изучить движение двух притягивающихся или отталкивающихся частиц, абсолютно упругий и неупругий центральные удары, абсолютно упругий и неупругий нецентральные удары, движение частицы в центрально-симметричном поле другой частицы, движение молекул газа, диффузия, движение планет вокруг Солнца, движение взаимодействующих частиц в однородном поле, движение взаимодействующих частиц в центрально-симметричном поле.

3. Алгоритм.

1. Задают число материальных точек , их массы , координаты и проекции начальных скоростей силовое поле

,

а также шаг по времени

2. Начало цикла по Дают приращение по времени: переменной присваивают значение

3. Определяют проекции равнодействующей всех внешних и внутренних сил, действующих на каждую ую материальную точку в момент , и записывают их в массивы.

4. В цикле переобозначают координаты всех материальных точек, записывая их в массивы

5. В цикле перебирают все материальные точки и определяют проекции ускорения, скорости и координаты для каждой из них в момент

По аналогичным формулам вычисляют проекции на ось OY. Результаты записывают в массивы

6. Стирают изображения материальных точек в предыдущий момент времени , координаты которых сохранены в массивах

7. На экране строят точки в следующий момент , либо рисуют графики или выводят результат в числовом виде.

8. Возвращение к операции 2. Если цикл по закончился, - выход из цикла.

4. Компьютерная программа. Ниже приведен код программы, которая моделирует движение 50 молекул газа в прямоугольном сосуде, находящемся в однородном гравитационном поле.

program PROGRAMMA4;

uses dos, crt, graph;

const N=50; dt=0.01;

var m,Fx,Fy,x,y,vx,vy,xx,yy : array[1..N] of real;

Gd, Gm, i, j : integer; ax, ay, F, l : real;

label Metka, metka1;

Procedure Init_Graph; {Инициализация графики}

begin

Gd:=Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');

if GraphResult <> grOk then Halt(1);

end;

Procedure Sila; {Вычисление действующих сил}

label Metka;

begin

For i:=1 to N do begin Fx[i]:=0; Fy[i]:=0; end;

For i:=1 to N do for j:=1 to N do begin

if j=i then goto Metka;

l:=sqrt(sqr(x[i]-x[j])+sqr(y[i]-y[j])); if l<2 then l:=2;

F:=-50000*m[i]*m[j]/sqr(l)+500000*m[i]*m[j]/sqr(l*l);

Fx[i]:=Fx[i]+F*(x[i]-x[j])/l;

Fy[i]:=Fy[i]+F*(y[i]-y[j])/l+m[i]*10;

Metka: end;

end;

Procedure Nach_uslov;

begin Randomize; {Задание случайных координат и скоростей}

for i:=1 to N do begin m[i]:=2;

x[i]:=random(280)+60; y[i]:=random(280)+60;

vy[i]:=random(30)-15; vx[i]:=random(30)-15;

end; end;

BEGIN

Init_Graph; Nach_uslov;

Repeat Sila;

for i:=1 to N do

begin

xx[i]:=x[i]; yy[i]:=y[i]; {Запись предыдущих координат}

ax:=Fx[i]/m[i]; ay:=Fy[i]/m[i]; {Вычисление ускорений,}

vx[i]:=vx[i]+ax*dt; vy[i]:=vy[i]+ay*dt; {скоростей,}

x[i]:=x[i]+vx[i]*dt; y[i]:=y[i]+vy[i]*dt; {координат}

if (x[i]<50)or(x[i]>350) then vx[i]:=-vx[i];{отражение}

if (y[i]<50)or(y[i]>350) then vy[i]:=-vy[i];{от стенок}

end; delay(500);

setcolor(8); for i:=1 to N do circle(round(xx[i]),round(yy[i]),2);

setcolor(15); for i:=1 to N do circle(round(x[i]),round(y[i]),2);

until KeyPressed; Repeat until keypressed; CloseGraph;

END.

5. Задания для студентов.

1. На основе компьютерной модели изучите движение двух или трех частиц, между которыми действуют силы притяжения.

2. Изучите движение двух материальных точек, между которыми действуют силы отталкивания. Промоделируйте центральный и нецентральный удар.

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

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

5. Промоделируйте движение нескольких планет и комет Солнечной системы, учитывая, что масса Солнца во много раз больше массы любой планеты.

6. Промоделируйте движение молекул газа в прямоугольном замкнутом сосуде. Учтите, что при соударении молекулы с горизонтальной (вертикальной) стенкой сосуда вертикальная (горизонтальная) проекция ее скорости меняет свой знак на противоположный.

7. Промоделируйте диффузию двух газов. Пусть вначале молекулы с массами заполняют левую половину сосуда, а молекулы с массами - правую. Задайте случайные значения скоростей молекул. Как изменяется концентрация молекул газов в сосуде с течением времени?

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

9. Промоделируйте движение молекул газа в гравитационном поле шара большой массы.

6. МОДЕЛИРОВАНИЕ КОЛЕБАНИЙ СВЯЗАННЫХ ОСЦИЛЛЯТОРОВ

1. Задача: Имеется система из осцилляторов с массами и жесткостью пружин , связанные между собой упругими связями с жесткостью . На каждый осциллятор действует сила вязкого трения Заданы начальные смещения и скорости всех осцилляторов. На отдельные осцилляторы действует вынуждающая сила

;

некоторые осцилляторы колеблются по закону

.

Известно, что крайние осцилляторы закреплены (свободны). Необходимо построить компьютерную модель колебаний связанных осцилляторов.

2. Теория. Рассмотрим цепочку осцилляторов с массами и коэффициентом жесткости , связанных между собой пружинами жесткостью . На каждый -ый осциллятор со стороны соседних -ого и -ого осцилляторов действует сила:

Кроме нее на осциллятор действует сила упругости и сила вязкого трения , где - скорость -ого осциллятора. По второму закону Ньютона, ускорение -ого осциллятора равно:

Определив ускорение , можно найти скорость и смещение -ой массы в следующий момент времени :

3. Алгоритм. Эта задача является частным случаем задачи о моделировании движения системы взаимодействующих между собой частиц.

1. Задают параметры системы и начальные условия: число осцилляторов , их массы , жесткость упругих связей , жесткость пружины внутри осциллятора , шаг по времени начальные смещение и скорости , вынуждающие силы, действующие на различные материальные точки.

2. Дают приращение по времени: переменной присваивают значение

3. В цикле перебирают все осцилляторов и для каждого из них определяют равнодействующую силу, ускорение, скорость и смещение в следующий момент

4. В цикле стирают предыдущее изображение каждого осциллятора и рисуют новое.

5. Возвращение к операции 2. Если цикл по закончился, - выход из цикла.

4. Компьютерная программа. Программа, представленная ниже, моделирует распространение импульса в случае, когда левый элемент среды совершает полколебания.

program PROGRAMMA5;

uses dos, crt, graph;

Const m=0.5; r=0.1; k=0.01;

dt=0.001; q=100; N=50;

Var teta, F, t, y : Real;

i,xx,vv,aa,FF,tt,Gd,Gm : Integer;

eta, ksi : array [0..N] of real;

Procedure Graph_Init;

begin

Gd:= Detect;

InitGraph(Gd, Gm, 'c:\bp\bgi');

if GraphResult <> grOk then Halt(1);

end;

Procedure Oscillator;

begin

F:=q*(ksi[i-1]-ksi[i])+q*(ksi[i+1]-ksi[i]);

teta:=(F-r*eta[i]-k*ksi[i])/m;

eta[i]:=eta[i]+teta*dt; ksi[i]:=ksi[i]+eta[i]*dt;

end;

BEGIN

Graph_Init;

Repeat begin t:=t+dt;

For i:=1 to N do begin

y:=ksi[i]; Oscillator;

if 5*t<3.142 then ksi[1]:=20*sin(5*t) else ksi[1]:=0;

ksi[N]:=0; {закреплен}

{ksi[N]:=ksi[N-1]; незакреплен}

setcolor(8); circle(10*i,240-round(y*10),2);

setcolor(15); circle(10*i,240-round(ksi[i]*10),2);

end; end;

until KeyPressed;

CloseGraph;

END.

5. Задания для студентов.

1. Промоделируйте колебания двух связанных осцилляторов. Рассмотрите случаи: 1) на один из них действует вынуждающая сила; 2) один из осцилляторов имеет начальное смещение; 3) один из осцилляторов имеет начальную скорость. Выполните компьютерные эксперименты при различных и .

2. Изучите колебания трех связанных осцилляторов, рассмотрев все перечисленные выше случаи, выполнив компьютерные эксперименты при различных и .

3. Промоделируйте колебания 50 осцилляторов, связанных упругими связями, в случае, когда на левый крайний осциллятор подействовала кратковременная сила. Рассмотрите случаи, когда правый крайний осциллятор закреплен и не закреплён.

6. Промоделируйте распространение импульса вдоль цепочки осцилляторов, связанных упругими связями, в случае, когда их масса или жесткость пружин, начиная с некоторого осциллятора, изменяется скачком. Изучите изменение фазы импульса при его отражении от "более плотной" и "менее плотной" среды.

7. Изучите распространение импульса и его отражение от открытого или закрытого конца струны (одномерной упругой среды), которая моделируется 50 связанными осцилляторами.

7. МОДЕЛИРОВАНИЕ РАСПРОСТРАНЕНИЯ ОДНОМЕРНОЙ ВОЛНЫ

1. Задача. Имеется непрерывная одномерная упругая среда. Заданы колебания отдельных элементов среды и скорость распространения возмущения. Необходимо рассчитать смещения элементов среды в последующие моменты времени.

2. Теория. Запишем волновое уравнение:

где - фазовая скорость волны. Скорость любой точки среды в момент составляет:

Разобьём среду на элементов равной длины смещение которых из положения равновесия обозначим через где - номер элемента. Скорость -го элемента в момент равна:

где - ускорение ого элемента:

Тогда смещение -ого элемента в момент времени из положения равновесия составляет

3. Алгоритм. Для построения "моментальной фотографии волны" в последующие моменты времени необходимо составить цикл по , в котором перебираются все элементы среды, вычисляются их смещения из положения равновесия. После этого стирается предыдущая моментальная фотография волны и строится новая. Этот цикл должен находиться внутри цикла по времени. Алгоритм модели представлен ниже.

1. Задают число элементов скорость распространения волны, а также уравнения колебаний отдельных элементов среды (например:

).

2. Начало цикла по Дают приращение по времени: переменной присваивают значение

3. Вычисляют смещение элементов среды, которые колеблются по заданному закону.

4. В цикле перебираются все элементы от 2 до . При этом записывают значения смещений в массив и вычисляют скорости элементов в момент времени

записывая их в массив

5. В цикле перебираются все элементы и вычисляются их смещения по формуле:

6. В цикле перебирают все элементы, стирают их предыдущие изображения и рисуют новые.

7. Возвращение к операции 2. Если цикл по закончился, -- выход из цикла.

4. Компьютерная программа. Предлагаемая программа моделирует прохождение и отражение импульса от "границы раздела двух сред".

program PROGRAMMA6;

uses crt, graph;

const n=200; h=1; dt=0.05;

var i, j, DriverVar,

ModeVar, ErrorCode : integer;

eta,xi,xxii : array[1..N] of real;

t, vv : real;

Procedure Graph_Init;

begin {- Инициализация графики -}

DriverVar:=Detect;

InitGraph(DriverVar,ModeVar,'c:\bp\bgi');

ErrorCode:=GraphResult;

if ErrorCode<>grOK then Halt(1);

end;

Procedure Raschet; {Расчет смещения}

begin for i:=2 to N-1 do

begin

xxii[i]:=xi[i];

if i<N/2 then vv:=4 else vv:=0.5;

eta[i]:=eta[i]+vv*(xi[i+1]-2*xi[i]+xi[i-1])/(h*h)*dt;

end;

for i:=2 to N-1 do xi[i]:=xi[i]+eta[i]*dt;

xi[N]:=0; {Конец закреплен}

{ xi[N]:=xi[N-1];}{ незакрепленный}

end;

Procedure Draw;

begin {- Вывод на экран -}

setcolor(black);

line(i*3-3,240-round(xxii[i-1]*50),i*3,240-round(xxii[i]*50));

setcolor(white);

line(i*3-3,240-round(xi[i-1]*50),i*3,240-round(xi[i]*50));

end;

BEGIN {- Основная программа -}

Graph_Init;

Repeat t:=t+dt;

if t<6.28 then xi[1]:=2*sin(t) else xi[1]:=0;

Raschet; For i:=1 to N do Draw;

until KeyPressed; CloseGraph;

END.

Рассмотренная выше компьютерная модель позволяет выполнить серию численных экспериментов и изучить следующие явления:

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

2) интерференция волн (одиночных импульсов, цугов), возникающая в результате отражения падающей волны либо излучения двух когерентных волн;

3) отражение и прохождение волны (одиночного импульса, цуга) через границу раздела двух сред;

4) изучение зависимости длины волны от частоты и скорости распространения;

5) наблюдение изменения фазы отраженной волны на при отражении от среды, в которой скорость волны меньше.

5. Задания для студентов.

1. Промоделируйте распространение волны и ее отражение от закрепленного (незакрепленного) правого конца среды в случае, когда ее левый элемент совершает гармонические колебания.

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

3. Пронаблюдайте суперпозицию волн, испускаемых двумя элементами, колеблющимися с равными (различными) частотами и отстоящими друг от друга на расстояние .

4. Промоделируйте возникновение стоячей волны при отражении гармонической волны от правого закрепленного (незакрепленного) конца шнура.

5. Изучите интерференцию двух цугов, распространяющихся навстречу.

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

7. Используя модель, изучите зависимость длины волны от частоты.



8. МОДЕЛИРОВАНИЕ ЯВЛЕНИЙ ПЕРЕНОСА

1. Задача: Имеется прямоугольная пластина, коэффициент температуропроводности которого зависит от координаты

.

Задано начальное распределение температуры

,

а также мощность и координаты источников тепла. Необходимо изучить изменение температуры различных точек пластины с течением времени.

2. Теория. Уравнение теплопроводности для пластины:

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

Приращение температуры равно:

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

Решая одномерную задачу теплопроводности по столбцам, получаем:

где - быстрота увеличения температуры элемента вследствие наличия в нем источника тепла, . Аналогичным образом решается уравнение теплопроводности по строкам:

Организовав два цикла по и , в которых перебираются все элементы пластины, вложенные в цикл по времени, можно рассчитать распределение температуры пластины в последующие моменты времени. Чтобы исключить влияние направления перебора элементов на результаты моделирования, следует организовать 4 цикла, в которых элементы пластины перебираются слева направо, справа налево, снизу вверх, сверху вниз.

3. Алгоритм.

1. Задают коэффициенты температуропроводности и вдоль осей координат и , начальное распределение температуры различных элементов пластины, координаты и мощности источников тепла . Считают, что

2. Путем последовательного перебора элементов -ого столбца от до с помощью конечно-разностного уравнения (1) пересчитывают их температуры в следующий дискретный момент времени. Повторяют расчеты для всех столбцов.

3. Повторяют ту же самую процедуру в противоположном направлении, перебирая элементы от до и используя уравнение (1).

4. Путем последовательного перебора элементов -ой строки от до с помощью конечно-разностного уравнения (2) пересчитывают их температуры в следующий дискретный момент времени. Повторяют расчеты для всех строк.

5. Повторяют ту же самую процедуру в противоположном направлении, перебирая элементы от до и используя уравнение (2).

6. Вводят текущее распределение температуры на экран, закрашивая элементы с различной температурой различными цветами.

7. Возвращение к операции 2. Если цикл по закончился, -- выход из цикла.

4. Компьютерная программа. Предлагаемая компьютерная программа позволяет исходя из начального распределения температуры и наличия источников тепла, рассчитать температуру различных точек прямоугольной пластины в дискретные моменты времени.

program PROGRAMMA7;

uses crt, graph;

const n=100; m=100; h=1; dt=0.2;

var ii,jj,kk,i,j,DriverVar, ModeVar, ErrorCode : integer;

t: array[1..N, 1..M] of real; q,a,b,bb :real; naprav,uslovie: boolean;

procedure Init; {---- Инициализация графики ----}

begin DriverVar:=Detect; InitGraph(DriverVar,ModeVar,'c:\bp\bgi');

ErrorCode:=GraphResult; if ErrorCode <> grOK then Halt(1); end;

procedure Param_sred; {---Коэффициент температуропроводности---}

begin if j<70 then a:=2 else a:=1; end;

procedure Istoch; {--- Источники тепла ---}

begin if ((i>45)and(i<55))and((j>70)and(j<75)) then q:=50 else q:=0; end;

procedure Nach_uslov; {-- Начальное распределение температуры --}

begin For i:=1 to N do For j:=1 to M do

begin uslovie:=((j<65)and(j>45)and(i>20)and(i<30))

or((j<45)and(j>35)and(i>50)and(i<60));

if uslovie=true then t[i,j]:=450 else t[i,j]:=1; end;

end;

procedure Raschet; {---- Расчет температуры ----}

begin Istoch; Param_sred;

t[i,j]:=t[i,j]+a*(t[i,j+1]-2*t[i,j]+t[i,j-1])*dt/(h*h)+q;

if naprav=true then t[i,j]:=

t[i,j]+a*(t[i+1,j]-2*t[i,j]+t[i-1,j])*dt/(h*h);

end;

procedure Draw;{---- Вывод на экран ----}

begin if t[i,j]>50 then setcolor(2);

if t[i,j]>300 then setcolor(12);

if (t[i,j]<300)and(t[i,j]>120) then setcolor(10);

if (t[i,j]<120)and(t[i,j]>70) then setcolor(3);

if (t[i,j]<70)and(t[i,j]>30) then setcolor(4);

if (t[i,j]<30)and(t[i,j]>20) then setcolor(5);

if (t[i,j]<20)and(t[i,j]>10) then setcolor(7);

if t[i,j]<10 then setcolor(15);

rectangle(i*5+50,j*5,i*5+54,j*5+4);

end;

BEGIN {---- Основная программа ----}

Init; Nach_uslov;

Repeat kk:=kk+1;

For i:=2 to N-1 do For j:=2 to M-1 do begin naprav:=false; Raschet; end;

For j:=2 to M-1 do For i:=2 to N-1 do begin naprav:= true; Raschet; end;

For i:=2 to N-1 do For jj:=2 to M-1 do

begin j:=M+1-jj; naprav:=true; Raschet; end;

For j:=2 to M-1 do For ii:=2 to N-1 do

begin i:=N+1-ii; naprav:=false; Raschet; end;

if kk/2=round(kk/2) then For i:=2 to N-1 do For j:=2 to M-1 do Draw;

until KeyPressed; CloseGraph;

END.

5. Задания для студентов.

1. Температура группы элементов, находящихся в центре стержня, достаточно высока. Постройте график зависимости температуры от координаты и исследуйте изменение распределения температуры вдоль стержня c течением времени, если коэффициент температуропроводности во всех точках одинаков.

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

3. Вблизи центра стержня имеется несколько источников тепла. Изучите изменение распределения температуры c течением времени, если стержень однороден.

4. Решите предыдущую задачу для случая, когда стержень неоднороден, то есть его коэффициент температуропроводности зависит от координаты.



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

6. Задайте источник тепла, мощность которого периодически изменяется с течением времени с очень низкой частотой. Промоделируйте тепловые волны.

7. Температура группы элементов, находящихся в центре пластины, достаточно высока. Исследуйте изменение распределения температуры c течением времени, если пластина однородная и изотропная.

8. Решите предыдущую задачу для случая, когда пластина неоднородна.

9. Промоделируйте нагревание неизотропной пластины источниками тепла, находящимися в центре.

10. Вблизи центра пластины имеется группа поглотителей тепла (источников тепла с отрицательной мощностью). Изучите изменение распределения температуры c течением времени.

11. Пластина с отверстием содержит источник тепла и поглотитель тепла. Изучите распределение температуры в различные моменты времени.

12. Температура группы элементов вблизи центра пластины поддерживается постоянной. Изучите распределение температуры, если пластина имеет источники тепла с положительной (отрицательной) мощностью.

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

9. МОДЕЛИРОВАНИЕ АВТОВОЛНОВЫХ ПРОЦЕССОВ

1. Задача: Имеется двумерная активная среда, состоящая из элементов, каждый из которых может находиться в трех различных состояниях: покое, возбуждении и рефрактерности. При отсутствии внешнего воздействия, элемент находится в состоянии покоя. В результате воздействия элемент переходит в возбужденное состояние, приобретая способность возбуждать соседние элементы. Через некоторое время после возбуждения элемент переключается в состояние рефрактерности, находясь в котором он не может быть возбужден. Затем элемент сам возвращается в исходное состояние покоя, то есть снова приобретает способность переходить в возбужденное состояние. Необходимо промоделировать процессы, происходящие в двумерной активной среде при различных параметрах среды и начальном распределении возбужденных элементов.

2. Теория. Рассмотрим обобщенную модель Винера-Розенблюта. Мысленно разобьём экран компьютера на элементы, определяемые индексами и образующими двумерную сеть. Пусть состояние каждого элемента описывается фазой , и концентрацией активатора , где - дискретный момент времени.

Если элемент находится в покое, то будем считать, что . Если вследствие близости возбужденных элементов концентрация активатора достигает порогового значения то элемент возбуждается и переходит в состояние 1. Затем на следующем шаге он переключается в состояние 2, затем - в состояние 3 и т.д., оставаясь при этом возбужденным. Достигнув состояния элемент переходит в состояние рефрактерности. Через () шагов после возбуждения элемент возвращается в состояние покоя. При этом реализуется следующее правило:

Если

Если

Если

Если

Будем считать, что при переходе из состояния в состояние покоя концентрация активатора становится равной . При наличии соседнего элемента, находящегося в возбужденном состоянии, она увеличивается на 1. Если ближайших соседей возбуждены, то на соответствующем шаге к предыдущему значению концентрации активатора прибавляется число возбужденных соседей:

Можно ограничиться учетом ближайших восьми соседних элементов.

3. Алгоритм. Для моделирования автоволновых процессов в активной среде необходимо составить цикл по времени, в котором вычисляются фазы элементов среды в последующие моменты времени и концентрация активатора, стирается предыдущее распределение возбужденных элементов и строится новое. Алгоритм модели представлен ниже.

1. Задают число элементов активной среды, ее параметры начальное распределение возбужденных элементов.

2. Начало цикла по Дают приращение по времени: переменной присваивают значение

3. Перебирают все элементы активной среды, определяя их фазы и концентрацию активатора в момент

4. Очищают экран и строят возбужденные элементы активной среды.

7. Возвращение к операции 2. Если цикл по закончился, -- выход из цикла.

4. Компьютерная программа. Ниже представлена программа, моделирующая активную среду и происходящие в ней процессы. В программе заданы начальные значения фазы всех элементов активной среды, а также имеется цикл по времени, в котором расчитываются значения в следующий момент и осуществляется графический вывод результатов на экран. Параметры среды , то есть каждый элемент кроме состояния покоя может находиться в 6 возбужденных состояниях и 7 состояниях рефрактерности. Пороговое значение концентрации активатора равно 5. Программа строит однорукавную волну, пейсмекер и препятствие.

Program PROGRAMMA8;

uses dos, crt, graph;

Const N=110; M=90; s=13; r=6; h=5;

Var y, yy, u : array [1..N,1..M] of integer;

ii, jj, j, k, Gd, Gm : integer; i : Longint;

Label met;

BEGIN

Gd:= Detect; InitGraph(Gd, Gm, 'c:\bp\bgi');

If GraphResult <> grOk then Halt(1);

setcolor(8); setbkcolor(15);

(* y[50,50]:=1; { Одиночная волна } *)

For j:=1 to 45 do { Однорукавная волна }

For i:=1 to 13 do y[40+i,j]:=i;

(* For j:=1 to M do { Двурукавная волна }

For i:=1 to 13 do begin y[40+i,j]:=i;

If j>40 then y[40+i,j]:=14-i; end; *)

Repeat

If k=round(k/20)*20 then y[30,30]:=1; {Пейсмекер 1}

(* If k=round(k/30)*30 then y[20,50]:=1; {Пейсмекер 2} *)

For i:=2 to N-1 do For j:=2 to M-1 do begin

If (y[i,j]>0) and (y[i,j]<s) then yy[i,j]:=y[i,j]+1;

If y[i,j]=s then begin yy[i,j]:=0; u[i,j]:=0; end;

If y[i,j] <> 0 then goto met;

For ii:=i-1 to i+1 do For jj:=j-1 to j+1 do begin

If (y[ii,jj]>0) and (y[ii,jj]<=r) then u[i,j]:=u[i,j]+1;

If u[i,j]>=h then yy[i,j]:=1; end;

met: end; Delay(2000); {Задержка}

cleardevice;

For i:=21 to 70 do begin

yy[i,60]:=0; yy[i,61]:=0; {Препятствие}

circle(6*i-10,500-6*60,3); circle(6*i-10,500-6*61,3); end;

For i:=1 to N do For j:=1 to M do

begin y[i,j]:=yy[i,j]; setcolor(12);

If (y[i,j]>=1) and (y[i,j]<=r) then circle(6*i-10,500-6*j,3);

setcolor(8);

If (y[i,j]>6) and (y[i,j]<=s) then circle(6*i-10,500-6*j,2);

end;

until KeyPressed;

CloseGraph;

END.

Размещено на Allbest.ru

...