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

При расчете физических процессов составляется математическая модель - система уравнений, описывающая зависимости между физическими величинами при некоторых упрощающих допущениях. Например, при движении точки вблизи поверхности Земли полагается ускорение свободного падения постоянным, не зависящим от высоты расположения точки над поверхностью. Для тел, движущихся с небольшой скоростью или в разряженной атмосфере, пренебрегают сопротивлением воздуха. Само точка часто заменяют материальной точкой, т. е. размерами точки пренебрегают. Физические процессы описываются, как правило системой дифференциальных уравнений, для решения которой применяют различные численные методы (модели). Широко используется метод конечных разностей, в котором бесконечно малые приращения переменных заменяют малыми (конечными) приращениями.

Например, изменение параметра времени представляют в виде:dt=t₂-t₁, а изменение функции "Х": dX(t) = X(t)-X(t-dt) = X(t₂)-X(t₁) = X₂-X₁.

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

Координаты точки X(t), Y(t) в некоторый момент времени "t" можно определить, зная координаты точки X(t-dt), Y(t-dt) в предыдущий момент времени "t-dt" и изменение (приращение) координат dX, dY:

X(t) = X(t-dt) + dX(t),

Y(t) = Y(t-dt) + dY(t).

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

dX(t) = Vx(t)dt,

dY(t) = Vy(t)dt.

Здесь Vx(t), Vy(t) - проекции скорости на оси координат.

Составляющие скорости Vx(t) и Vy(t) можно вычислить по формулам:

Vx(t) = Vx(t-dt) + Ax(t)*dt,

Vy(t) = Vy(t-dt) + Ay(t)*dt.

Здесь Ax(t), Ay(t) - проекции ускорения на оси координат.

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

Приведем алгоритм расчета траектории движения точки:

1. Определяем силы, действующие на точка, и находим проекции ускорения на оси координат.В общем случае ускорение точки зависит от многих факторов и в момент времени t задается как функция от времени, скорости и координат точки:

Ax:= Fx(Vx, Vy, X, Y, t); Ay:= Fy(Vx, Vy, X, Y, t);

Где Vx, Vy, Ax, Ay - проекции скорости и ускорения.

2.Задаем начальное положение точки - координаты X[1], Y[1] и начальную скорость и ускорение в виде проекций на оси координат:

X[1]:= X0; Y[1]:= Y0; Vx[1]:= V*cos(fi); Vy[1]:= V*sin(fi);

Ax[1]:= Fx(Vx[1], Vy[1], X[1], Y[1], t[1]);

Ay[1]:= Fy(Vx[1], Vy[1], X[1], Y[1], t[1]);

Где V - начальная скорость точки, fi - угол наклона вектора скорости к оси Х.

3. Задаем временной шагdt и разбиваем весь временной интервал на N участков. При равномерной разбивке приращение времени определяется по формуле:

dt:= (t[N]-t[1])/(N-1); Здесь (t[N] - t[1]) - время движения точки.

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

4.Вычисляем массивы скорости, ускорения и координат точки:

For i:= 2 to N do begin

Vx[i]:= Vx[i-1] + Ax[i-1]*dt;

Vy[i]:= Vy[i-1] + Ay[i-1]*dt;

X[i]:= X[i-1] + 0.5*(Vx[i-1] + Vx[i])*dt;

Y[i]:= Y[i-1] + 0.5*(Vy[i-1] + Vy[i])*dt;

Ax[i]:= Fx(Vx[i], Vy[i], X[i], Y[i], t[i]);

Ay[i]:= Fy(Vx[i], Vy[i], X[i], Y[i], t[i]);

{ уточняем скорость точки в расчетной точке }

VX[i]:= VX[i-1] + 0.5*(Ax[i-1] + Ax[i])*dt;

VY[i]:= VY[i-1] + 0.5*(Ay[i-1] + Ay[i])*dt;

end;

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

5. Строим траекторию движения точки. Здесь удобно использовать процедуры из библиотеки построения графиков GR_F. Следует определить расчетную область и область рисования траектории на экране. Траектория на экране рисуется процедурой: PutPixel_G(X[i], Y[i], N);

Для тестирования работы алгоритма рассмотрим задачу расчета траектории точки, движущегося из точки с координатами X, Y с начальной скоростью Vx, Vy под действием сил, вызывающих ускорение точки Ax, Ay. Следуя пунктам 1. . 5 приведенного выше алгоритма необходимо рассчитать траекторию движения точки и сравнить с траекторией точки, описанной аналитической зависимостью X(t), Y(t).

Практическое задание N 2. 22

Рассчитать разностным моделированием и по аналитической зависимости траектории точки. Параметр a = 10, b = 5, 0<= t <=4*pi, N=500. Построить траектории точки.

N X₁ Y₁ Vx₁ Vy₁ Axi Ayi X(t) Y(t)

1 0 0 0 b 2*a -y a*t² b*sin(t)

2 0 0 a b 0 -y a*t b*sin(t)

3 1 0 1 1 -2*y 2*x e^t* cos(t) e^t*sin(t)

4 a 0 0 0 -x x*b/a a* cos(t) b*(1-cos(t))

5 a b 0 0 -4*x y a* cos(2*t) b*cos(t)

6 0 0 0 b 2*a 0 a*t² b*t

7 2*a 0 0 a x 0 a*(e^t + e^-t) a*t

8 0 b a 0 -x -y a* sin(t) b*cos(t)

Y V F, * V0 g fi 0 X

Рассмотрим задачу расчета траектории снаряда, движущегося с начальной скоростью "V₀" под углом "fi" к горизонту с учетом сил сопротивления воздуха, пропорциональных скорости снаряда. Проекции ускорений определим в виде функций:

FUNCTION Fx(Vx, kc: real): real; begin Fx:= - kc*Vx end;

FUNCTION Fy(Vy, kc: real): real; begin Fy:= - kc*Vy - g end;

Где kc - коэффициент сопротивления воздуха,

g = 9. 81, м/с - ускорение свободного падения у поверхности Земли.

Поскольку время подлета снаряда к цели неизвестно, то параметр "dt" выбирается приближенно, например, исходя из максимального времени полета снаряда над горизонтальной поверхностью без учета сопротивления воздуха: tмах= 2*V*sin(fi)/g. Для N = 500, dt = t/500. При решении конкретных задач процесс расчета прекращается при достижении снарядом цели, либо при ограничениях по статическим координатам, например:

i:= 1;

REPEAT i:=i+1;

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

Until ( cc = GetPixel_G(X[i], Y[i])) or ( Y[i] < 0 ) or ( i = N );

Здесь cc - цвет пикселов цели, Y[i] < 0 - ограничение по горизонтальной поверхности, i = N - ограничение по размеру массива. В случае преждевременного завершения полета снаряда необходимо увеличить dt или параметр N.

Практическое задание N 2. 23

1. Рассчитать разностным моделированием и по аналитической зависимости траектории полета снаряда без учета сопротивления воздуха. Построить траектории полета снаряда. Начальная скорость V₀=1000, м/с, угол fi=450. Аналитическая зависимость имеет вид:

X = V₀*t*cos(fi); Y = V₀*t*sin(fi) - g*t²/2;

2. Рассчитать разностным моделированием и по аналитической зависимости траектории полета снаряда с учетом сопротивления воздуха, пропорциональным скорости снаряда. Построить траектории полета снаряда. Начальная скорость V₀=3000, м/с, угол fi = 45⁰. Коэффициент сопротивления воздуха kc = 0. 01,с^-1.

Аналитическая зависимость имеет вид:

X=V₀*cos(fi)*(1-e^(-kc*t))/kc; Y=(V₀*sin(fi)+g/kc)*(1-e^(-kc*t))/kc-g*t/kc;

3. Рассчитать разностным моделированием траектории полета снаряда с учетом сопротивления воздуха, пропорциональным квадрату скорости снаряда. Коэффициент сопротивления воздуха kc₁ = kc² . Построить совместно траектории полета снаряда для п. 1, 2, 3. Начальная скорость V₀ = 3000, м/с, угол fi = 45⁰.

4. Составить программу поражения неподвижной цели при kc₁ = kc². Изменяя в цикле угол fi на небольшую величину, определить в программе угол при котором будет поражена цель - небольшой прямоугольник с координатами вершин (x1, y1) и (x2, y2). Построить все траектории полета снаряда.

Примечание к п. 1. . 4: Выводить на экран исходные данные: V₀, fi, kc, а также наибольшую высоту и дальность полета снаряда.

Рассмотрим задачу расчета траектории космического тела, в поле тяготения планеты без учета сил сопротивления. В начальный момент времени тело движется на высоте "Н" со скоростью "V₀", направленной по касательной к окружности радиуса R₀. Поскольку движение спутника вокруг планеты достаточно продолжительно, то не целесообразно запоминать в оперативной памяти все параметры (координаты, скорости и ускорения) в каждый момент времени. Обычно эти параметры, записываются в файл на диск при вычислениях через некоторые моменты времени, а траекторию строят сразу, либо запуская отдельную программу, считывающую данные из файла. Расчетная область задается исходя из оценочных расчетов. Для спутника, движущегося вокруг Земли, можно принять:

Xmin= Ymin= -Kv*R₀, Xmax= Ymax= Kv*R₀,

Здесь R₀ = (Rz+H), Rz=6. 37*10⁶, м. - радиус Земли.

Kv=1. 5 при V₀ <= W₁; Kv=10 при W₁ < V₀ < W₂; Kv=20 при V >= V₂.

W₁ = Rz*Ö(g/R₀) - первая космическая скорость,

W₂ = Ö2* W₁ - вторая космическая скорость.

Параметр "dt" можно определить приближенно по формуле: dt=T/N,

где T= 6. 28*Rz/W₁ - время оборота спутника вокруг Земли, N=300.

Расстояние от спутника до центра планеты определяется через координаты:

function R(x, y: double): double; begin R:= sqrt(x*x + y*y) end;

Проекции ускорений определим в виде функции:

function FA(x,r,kz: double):double; begin FA:= -kz*x/(r*r*r) end;

Здесь kz = 4. E+14 для Земли (в системе СИ).

Пусть в начальный момент времени известны координаты спутника:

x₁ = R₀; y₁ = 0; r₁ = R(x₁, y₁);

скорость: Vx₁ = 0; Vy₁ = V₀;

и ускорение: Ax₁ = FA(X₁, r₁, kz); Ay₁ = FA(Y₁, r₁, kz);

Отметим, что скорость в начальный момент времени направлена по касательной к окружности радиуса r₁.

Для записи алгоритма расчета траектории необходимо знание параметров в двух соседних точках, например, в точке "1" - для предшествующего момента времени и в точке "2" - для расчетного момента времени. Расчет производим в цикле с одновременным выводом траектории движения спутника на экран до тех пор пока выполняется ограничение по радиусу траектории или не нажата любая клавиша.

While ( r1< Xmax ) or ( r1> Rz ) or ( not keyPressed ) do begin

Vx2:= Vx1 + Ax1*dt; Vy2:= Vy1 + Ay1*dt;

X2:= X1 + 0.5*(Vx1 + Vx2)*dt;

Y2:= Y1 + 0.5*(Vy1 + Vy2)*dt; r2:= R(x2, y2);

Ax2:=FA(X2, r2, kz);

Ay2:=FA(Y2, r2, kz);

Vx2:= Vx1 + 0.5*(Ax1 + Ax2)*dt;{ уточняем скорость }

Vy2:= Vy1 + 0.5*(Ay1 + Ay2)*dt;

{ Переопределяем значения параметров в точке }

x1:= x2; y1:= y2; r1:= r2;

Vx1:= Vx2; Vy1:= Vy2; Ax1:= Ax2; Ay1:= Ay2

PutPixel_G(x1,y1,c); { Строим траекторию движения точки, c - цвет точки }

end;

Практическое задание N 2. 24

1. Рассчитать разностным моделированием и по аналитической зависимости траектории полета спутника Земли. Аналитическая зависимость имеет вид:

r = P/(1 + e*cos(fi));

где e = P/R₀ - 1; P = (V₀* R₀/Rz)²/g ; 0 <= fi = 2*Pi.

В начальный момент времени известны координаты спутника: x₁ = R₀; y₁ = 0;

и скорость: Vx₁ = 0; Vy₁ = V₀; Рассмотреть случаи:

1_1. Начальная скорость V₀ <= W₁, высота H = 300000, м.

1_2. Начальная скорость W₁ <= V₀ < W₂, высота H = 400000, м.

1_2. Начальная скорость V₀ >= W₂, высота H = 500000, м.

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

2. Рассчитать разностным моделированием и построить траектории полета спутника вокруг двух планет (типа “Земля”), при V₀ < W₂, в случаях:

1) V₀ Rz Rz 2) Rz V₀ Rz

3. Рассчитать разностным моделированием и построить траектории полета двух планет типа “Земля” и их центра масс, при V₀ < W₂, в случаях:

V₀ V₀ V₀

1) 20 *Rz 2) 20 *Rz

V₀

Рассмотрим задачу расчета траектории точки переменной массы, движущегося под действием реактивной тяги. Движение точки в этом случае описывается уравнением Мещерского:

A = (U/M)*(dM/dt) + F/M

Где A - ускорение точки, M - масса точки.

U - скорость реактивной струи относительно точки,

F - результирующая внешних сил, действующих на точку,

Учитывая, что F = kz*M/r²- сила притяжения направлена к центру Земли, а P = U*(dM/dt)- реактивная сила двигателя (тяга) направлена по касательной к траектории движения, определяем проекции ускорения на оси координат:

Ax = P*Vx/(M*V) - kz*x/(r³); Ay = P*Vy/(M*V) - kz*y/(r³);

Где V = Ö(Vx² + Vy² )- скорость точки,

r = Ö( x² + y² ) - расстояние до центра Земли,

Vx , Vy - проекции скорости точки на оси координат, x, y - координаты точки.

Полагая расход топливаz = dM/dtпостоянным, массу точки можно определить по формуле: M = M₀ - z*t; при t < Tk,

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

Практическое задание N 2. 25

1. Построить десять траекторий полета баллистической ракеты, рассчитанных разностным моделированием. Начальная скорость V₀=1,м/с, тяга двигателя P=2. 5Е6,н, стартовая масса M₀ = 1. 5Е5, кг, расход топлива z= 700, кг/с, время работы двигателя Tk = 200, с.

2. Построить траектории полета двухступенчатой баллистической ракеты, рассчитанные разностным моделированием. Начальная скорость V₀ = 1,м/с, стартовая масса M₀ = 3Е5, кг, для первой ступени: тяга P₁ =5Е6, н, расход топлива z₁= 1700, кг/с, время работы двигателя Tk₁ = 130, с. Для второй ступени: тяга P₂ = 1. 1Е6, н, расход топлива z₂= 300, кг/с, время работы двигателя Tk₂ = 230, с.

Примечание к п. 1, 2: сопротивление воздуха и вращение Земли не учитывать. Угол запуска ракеты к горизонту = 90⁰ -N*0. 002⁰, где N= 1, 2, 3, ..., 10. Во время работы двигателя dt=0. 05, c, затем dt=0. 5, c.

3. Построить траекторию полета спутника Земли при включении двигателя, рассчитанную разностным моделированием. Начальные условия на высоте H=400000 м принять следующие: скорость V₀=W₁ и направлена по касательной к окружности, M₀=11000, кг, тяга двигателя P=4Е5, н, расход топлива z=100, кг/с, время работы двигателя Tk = 70, с. Рассчитать скорость спутника при работе двигателя по формуле Циолковского: V = V₀ + U*ln(M₀/M), где U = P/z.

Через каждые 10 секунд выводить на экран время полета спутника и скорость.

Рассмотрим задачу расчета траектории точки, прикрепленной к упругой нити, и движущейся с начальной скоростью "V₁" под углом "fi" к оси "x" из точки с координатами (x₁, y₁), без учета сил сопротивления воздуха. Эта задача моделирует известную игрушку - мяч, привязанный на резинке.

Пусть точка имеет массу "M", длина нити "L". Полагаем, что нить невесома и абсолютно упруга. Коэффициент упругости "Kn".

Оси координат проведем через точку закрепления нити вверх и влево. Расчетную область ограничим: X_min = Y_min = -Lm, X_max = Y_max = Lm,

где Lm = abs(V₁* Ö(M/Kn)) + Ö(x12 + y12) + L + 2*M*g/Kn.

Y V1   x,y     0 X

Период свободных колебаний груза,

подвешенного на упругой нити:

T = 6, 28* Ö(M/Kn). Примем dt = T/300.

Проекции ускорения определяются как дискретная функция расстояния " r " от начала координат до точки закрепления нити: если r <= L, то ускорение от сил упругости равно нулю, в остальных случаях:

Ax = -x*Ky*dr/(r*M);

Ay = -y*Ky*dr/(r*M) - 9.81; где dr = (r-L) > 0.

Проекцию ускорения на ось “Х” от сил упругости, запишем в виде функции:

FUNCTION FA(x, r, L, Kn, M: double): double;

begin if (r-L)>0 then FA:= -x*Kn*(r-L)/(r*M) else FA:= 0 end;

Аналогичная функция составляется для проекции ускорения на ось “У”. Методика расчета соответствует приведенной для движения спутника в поле тяготения планеты.

Практическое задание N 2. 26

1. Построить траекторию движения мяча, подвешенного на упругой нити в вязкой среде, рассчитанную разностным моделированием. Сопротивление среды пропорционально скорости движения мяча: kc=0. 01, с^-1. Нить закреплена в центре квадрата со стороной 2*Lm, длина нити L=1, м, коэффициент упругости Kn=5, н/м. Масса мяча M=0. 2, кг. Мяч начинает движение из точки с координатами x₁=-0. 5*L, y₁=0, со скоростью V₁=10, м/с, под углом 45⁰.

2. Построить траекторию движения мяча, подвешенного на упругой нити в квадратной коробке, рассчитанную разностным моделированием, с учетом уменьшения нормальной составляющей скорости на 20% при отражении мяча от стенки. Сопротивление среды пропорционально скорости движения мяча: kc=0. 05, с^-1. Нить длиной L=1, м, закреплена в центре квадрата со стороной a=1. 5*L. Коэффициент упругости Kn=5, н/м, масса мяча M=0. 1, кг. Мяч начинает движение из точки с координатами x₁=-L, y₁=0, со скоростью V₁x=1, м/с, V₁y=5, м/с.

2. 4. Моделирование многовариантных задач с использованием графов

А   1 4 2 3 В

Рассмотрим "классический" пример многовариантной задачи. Пусть пункты A и B связаны между собой дорогами, могущими проходить также через пункты 1, 2, 3,..., N. В общем случае каждый пункт связан дорогами со всеми остальными. В частном случае некоторые связи (дороги) отсутствуют. Схематически эти пункты и связи можно изобразить в виде графа.

Графом называется совокупность узлов (пункты A, B, 1, 2, . . . , N) и связывающих их ребер (дорог). Маршрутом движения называется последовательность связанных ребрами узлов. В дальнейшем будем рассматривать те маршруты движения, которые всегда начинаются из пункта A и заканчиваются в пункте B. Причем пункты A и B на маршруте повторяться не могут. Например : А-1-4-В.

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

Пусть узел A имеет номер "0", а узел B - номер "N+1". Рассмотрим общий случай: каждый пункт связан со всеми остальными. Обозначим M - число промежуточных узлов на маршруте.

При М = 0 маршрут может проходить только из узла "0" в узел "N+ 1".

При М = 1 маршрут проходит через один из узлов: j1= 1, либо j1= 2, .., либо j1= N.

При М = 2 маршрут проходит через два узла, причем первый из них может иметь номер: j1=1, либо j1=2, ... либо j1=N, а второй - номер: j2=1, либо j2=2, ... либо j2=N, т. е. возможно N² маршрутов. Графически все маршруты можно представить в виде:

A M=1 A M=2

1 . . . j1 . . . N

1 2 3 ... j1 ... N 1 2 3 ... j2 . N 1 2 3 ... j2 ... N 1 2 3 ... j2 .. N

B B

Таким образом, число маршрутов равно N^M и время перебора маршрутов при больших значениях N и M очень быстро растет.

При постановке задачи нахождения маршрутов указывается значение M - наименьшее число узлов на маршруте, M1 - наибольшее число узлов на маршруте. Причем 1<=M<=M1. Например, пусть на графе имеется три узла N=3 и необходимо составить маршруты, проходящие через два узла, т. е. M=2, M1=2. Тогда в общем случае имеются маршруты:

A

0-1-1-4; 0-2-1-4; 0-3-1-4; односторонняя связь

0-1-2-4; 0-2-2-4; 0-3-2-4; 1 2 3

0-1-3-4; 0-2-3-4; 0-3-3-4; двусторонняя связь

B

Постановка задачи нахождения маршрутов включает определение матрицы коэффициентов aij, характеризующих связи между узлами i и j. Связь узла A задается коэффициентами a₀j, узла В - коэффициентами ai_{N+ 1}. Матрица имеет вид:

a₁₁ a₁₂ a₁₃ ... a_1N Если aij = aji = 0, то связь

a₂₁ a₂₂ a₂₃ ... a_2N между узлами i и j отсутствует.

a₃₁ a₃₂ a₃₃ ... a_3N Если aij=0 и aji<>0, то связь

........................... .между узлами i и j односторонняя.

a_N1 a_N2 a_N3 ... a_NN Если aij<>0 и aji<>0, то связь

между узлами i и j двусторонняя.

Если aij = aji при i =1, 2, . . , N; j = 1, 2, . . , N, то матрица симметричная.

Если aij = 0 при j =1, 2, . . , N; i > j, то матрица треугольная.

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

Введем линейный массив "Y", коэффициенты которого обозначают номера узлов графа через которые проходит маршрут, а индексы показывают номер пункта по порядку следования на маршруте. Операторы по перебору маршрутов имеют вид:

Y[0]:=0;{номер узла "А" графа}

repeat{цикл по числу узлов на маршруте}

for j:= 1 to M do Y[j]:=1;{начальные номера узлов на маршруте}

Y[M+1]:=N+1;{номер узла "B" графа}

repeat{цикл по перебору номеров узлов на маршруте}

for j:=1 to M+1 do if a[y[j-1],y[j]]=0 then goto METKA;{проверка}

{связей}

{****** здесь ставятся операторы фильтра ************}

{****** . .. . . . . . . . . . . . . . . . . . . . . . . . . . . . . ************}

for j:=0 to M+1 do write('-', Y[j]); writeln;{вывод маршрута}

METKA: Y[1]:=Y[1]+1; {изменение номера узла первого пункта на маршруте}

for j:=1 to M-1 do{определяем номера узлов на маршруте}

if Y[j]>N then begin Y[j]:=1; Y[j+1]:=Y[j+1]+1 end else Break;

until Y[M]=N+1;

M:=M+1

until M>M1;

В начале программы задается возможный маршрут 0-1-1-1-. . . -1-N+1 для заданного значения M>0. Проверяется наличие связей и ставятся фильтры для определения маршрута. Затем увеличивается номер узла первого пункта по порядку следования на маршруте: 0-2-1-1-. . . -1-N+1 и т. д. до 0-N-1-1-. . . -1-N+1. При превышении номера узла значения N, номер узла сбрасывается до единицы, а номер следующего узла увеличивается на единицу: 0-1-2-1-. . . -1-N+1 и снова увеличивается номер узла первого пункта до значения N: 0-N-2-1-. . . -1-N+1 и далее сбрасывается до единицы с увеличением номера следующего узла: 0-1-3-1-. . . -1-N+1. После (N-1)-го сброса и увеличения значения узла первого пункта до N получим маршрут: 0-N-N-1-. . . -1-N+1 и далее: 0-1-1-2-. . . -1-N+1. Таким образом, происходит перебор всех возможных маршрутов до 0-N-N-N-. . . -N-N+1. После этого рассматриваются маршруты для M=M+1 включая M=M1. Отметим, что при необходимости маршрут 0-N+1 для M=0 нужно рассмотреть отдельно.

При решении конкретных задач необходимо определить значение коэффициентов aij матрицы связи и установить необходимые фильтры.

Рассмотрим задачуопределения стоимости маршрутов из A в B.

1. ) Зададим стоимость проезда из узла i в узел j:

for i:=0 to N+1 do for j:=i to N+1 do a[i,j]:=Random(X); {X-дано}

for i:=0 to N+1 do a[i,i]:=0;{ движение внутри узла запрещено}

for i:=0 to N+1 do for j:=i to N+1 do a[j,i]:=a[i,j]; {связи }

{двусторонние и равнозначные}

2). Матрицу связей можно вывести на экран для проверки. При выводе маршрута на экран или в файл можно выводить также значение стоимости маршрута.

S:=0; for m:=1 to M1+1 do S:=S+a[y[m-1],y[m]]; {стоимость маршрута}

1 2 3   4 5 6   7 8 9

Рассмотрим задачу расстановки мин на прямоугольном поле размером Nx*Ny. При этом M=M1=N=Nx*Ny и все узлы должны быть пройдены без повторений. Расстановка начинается из узла с заданным номером NH и может закончиться в узлах на верхней границе.

1) Определим матрицу связей:

for i:=0 to N+1 do for j:=1 to N+1 do a[i,j]:=0;

for i:=1 to N-1 do begin a[i,i+1]:=1; a[i+1,i]:=1 end;{связи по гориз}

for j:=1 to Ny-1 do begin k:=Nx*j; a[k,k+1]:=0; a[k+1,k]:=0 end;

for i:=1 to Nx do for j:=1 to Ny-1 do{связи по вертикали}

begin k:=Nx*(j-1)+i; a[k,k+Nx]:=1; a[k+Nx,k]:=1 end;

a[0,NH]:=2;{ NH - узел связи c узлом 0}

for i:=1 to Nx do a[i,N+1]:=3;{ 1, . . , Nx - узлы связи c узлом N+1}

2). Установим фильтр, запрещающий возврат в узел на маршруте:

for k:=1 to M do c[y[k]]:=0; for k:=1 to M do

begin c[y[k]]:=c[y[k]]+1; if c[y[k]]=1 then goto METKA end;

Здесь производится суммирование повторяющихся номеров узлов на маршруте. При совпадении номера узла значение счетчика c[y[k]]=1 -маршрут не рассматривается.

Рассмотрим задачу загрузки N - видов коробокв машину. Задается число коробок каждого вида: Ki, их вес Mi и объем Vi, где i=1, 2, . . , N. Ограничения могут быть по общему весу и объему. Число узлов графа равно N. Число узлов на маршруте M=1, М1=K₁+K₂+. . . +K_N. Интервал М-М1 можно уменьшить просчитав наибольшее допустимое по весу и объему число коробок KMi каждого вида загружаемых в машину (KMi<=Ki). Тогда М = min(KMi), а М1 = max(KMi). Поскольку порядок загрузки не имеет значения, то все связи односторонние. 0

1 2 ... k ... N N+1

1) Определим матрицу связей:

for i:=0 to N+1 do for j:=i to N+1 do a[j,i]:=0; {нижний треугольник}

for i:=0 to N+1 do for j:=i to N+1 do a[i,j]:=1;{верхний треугольник}

2) Определение числа коробок каждого вида аналогично суммированию повторяющихся номеров узлов на маршруте.

Практическое задание N 2. 27

1) Вывести в файл стоимость маршрутов без повторяющихся узлов при N=4, M=3, M1=4, Х=9. Определить номера маршрутов с наименьшей и наибольшей стоимостью

для разных значений М.

2) Вывести символами псевдографики в текстовом режиме маршруты движения в прямоугольнике 2х4, либо 4х2. Начало движения при NH=8.

3) Вывести общий вес и число коробок каждого из 3-х видов, загружаемых в машину. Задать веса функцией Random(50)+50; Установить фильтр по общему весу G<900. Общее число коробок: M=10, M1=12.

Не нашли, что искали? Воспользуйтесь поиском по сайту: