Будко В.Н. Сборник описаний задач специального лабораторного практикума Компьютерный эксперимент темы 3, 4, 5
Подождите немного. Документ загружается.
Рассмотрим динамику измененния чисенности популяции (на единице площади) в процессе
рождения и гибели индивидов популяции при взаимодействии с окружающей средой.
Обозначим численность (или плотность) популяции в момент t : x(t). Относительный прирост
Мальтус предложил простую гипотезу о пропорциональности относительного прироста величине
популяции :
час,минута , и т.п./получим дискретную модель популяции по Мальтусу : x i+1=(k+1)x i.
Ферхлюст и Пирл предложили уточнить модель : определять коэф . прироста как разницу
коэф . рождаемости p>=0 и смертности c<=1 :k=p-c и считать, что коэф . смертности линейно зависит
от величины популяции :
x xпор x
Функция c(x) для этой модели показана на рисунке .
Представляет интерес модель, в которой коэф . смертности начинает увеличиваться лишь после
того как плотность популяции станет больше некоторого порога .
Программа моделирования : это может быть, например :
--- поиск стационарного состояния популяции ,
--- поиск момента , когда популяция вырождается, вымирает,
--- определение значений “a” и ”b” при которых наступают периодические колебания численности
популяции,
--- и другие задачи , см.[11].
5.2.Модель хищник-жертва.
Имеются популяции двух видов, один из которых питается другим. Поведение этой системы
описывается соотношениями :
где : xi –численность (плотность) жертв ,
yi – хищников,
a
1,
b
1
(
a
2,
b
2
)
коэф . рождаемости и смертности жертв (хищников),
t
txttx
x
будет
t
времени
отрезок
за
x
∆
−∆+
=
∆
)()(
ˆ
)(
:
)
(
δ
δ
)(
)()(
.
.
0
)
(
)
(
txk
t
txttx
прироста
коэф
k
t
x
k
x
⋅=
∆
−∆+
−
≥
⋅
=
δ
,/1 суткиtеденицеусловнойравнымвременициидискретизаtотрезокПолагая
=
∆
∆
xbc
⋅
+
=
α
2
1
2
1
)1(
,)()()(;
ˆ
,
ˆ
;
iii
iiii
bxxax
bxaxxxxbatxttxьноСледовател
p
a
x
b
a
x
b
p
k
Тогда
−⋅+=
−+=⋅−=−∆+
−
=
⋅
−
=
⋅
−
−
=
+
+
α
α
с
с
α
α
>−+
<
=
⋅ порпор
пор
xxbxx
xx
c
,(
,
)α
α
iiiiii
iiiiii
yxqybyayy
cycxyxqxbxaxx
2
2
221
2010,1
2
111 ,
−++=
==−+=
+
−+
q1- коэф . защиты жертв ,
q2 – коэф . прожорливости хищников.
Программа моделирования :
--- исследовать графики динамики каждого вида.
--- выяснить характер зависимости периодических колебаний численности жертв и хищников от
a,b,c,q.
--- и другие задачи , см.[11].
5.3.Динамика населения региона с учётом распределения по возрастным группам.
Пкажем как строится такая и им подобные модели .
Обозначим плотность (по площади) или число жителей региона x(t) .Обозначим среднюю
продолжительность жизни человека T лет (целое число ). Например, T=80 лет. Разобьём всё
население на “T” возрастных групп с интервалом в 1 год жизни . Тогда количество жителей ,
возраст которых в момент анализа равен t, t=1,2,3,… ,T составит массив N(t).Выделим из него три
подмножества. Подмножество детей :
Ng=N(i) ,i=1,2,3,… ,I
где I максимальный возраст этой групы (например, I=17 лет).Подмножество родителей :
Np=N(j), j=I+1,I+2,… ,J
где J – максимальный возраст в группе родителей (например,J=40 лет).
Подмножество пожилых :Nп=N(k) , k=J+!,J+2,… ,T.
Понимаем , что возрастная группа N(T) вся умирает через год естественной биологической
смертью.
Для подмножества родителей вводим характеристику : каждая семья за интервал лет от I+1 до J
имеет D детей.Следовательно для математичской модели можно считать, что в среднем каждый
родитель за интервал лет своей жизни от I+1 до J будет иметь D/2 детей.
Далее, для каждого подмножества вводим свои коэффициенты смертности (от болезней,
катастроф и других причин) :Cg,Cp,Cп<1.
Алгоритм модели .
1.Ввести исходные данные дали модели :T,I,J,D,Cg,Cp,Cп .
2.Ввести массив исходного распределения населения по возрастам N(T).Например, равномерное
распределение : N(T)=N0. Константу N0 удобно взять N0=1, понимая под этой условной
единицей,скажем, 1000 или 100000 жителей .
3.Вывести на экран график исходного распределения числа жителей по возрастам для начального
число жителей.(М g,Мр,Мп,М=Мg+Мр+Мп).
N М =М g+Мр+Мп
1 2 … .I I+1… . J J+1… .T t
Анализ динамики изменения населения по возрастам проводим в текущем дискретном времени
с интервалом в 1 год.
За каждый текущей год происходят следующие события.
1/ Так как каждый родитель (в среднем) за интервал жизни в [J – (I+1)] лет будет иметь D/2
детей,то , опять же в среднем , для динамики математической модели можно считать, что за 1 год в
каждой возрастной группе подмножества родителей появилось детей :
общееивподмножестжителейколичестваначальногоивременимоента 0
=
τ
0
=
τ
∑
=
=
I
t
tNМ g
1
)(
∑
+=
=
J
It
tNМр
1
)(
∑
+=
=
T
Jt
tNМп
1
)(
τ
Сумму детей появившихся за этот год у всех родителей затем следует присвоить переменной N(1).
2/ Количество жителей каждой возрастной группы уменьшилось на число умерших в этой
возрастной группе . Поэтому надо всем переменным N(t), t=1,2,… ,(T-1) присвоить новые значения :
N(t):=N(t) – CN(t)
со своими Cg,Cp,Cп для каждого подмножества возрастов .
3/ За этот же год анализа все возрастные группы постарели на 1 год. Поэтому надо все
переменные N(t) “сдвинуть на 1 год” , т.е . произвести переприсваивание :
N(k):=N(k – 1) , k=T,T-1,T-2,… ,4,3,2
При этом старшая возрастная группа N(t) после переприсваивания автоматически удаляется из
нового ряда N(t) – естественная биологическая смерть.
4/ Итак, следующими пунктами алгоритма должны быть циклы , реализующие происходящие
каждый год события.
5/ Вывод на экран нового распределения числа жителей по возрастам и вывод текущих значений
6/ Пауза 2-5 секунд для наблюдения графика.
Программа должна предусматривать возможность восприятия с клавиатуры во время паузы
команды останова анализа и приёма новых исходных данных для продолжения анализа динамики
населения в изменяющихся условиях .
7/ Возврат к пункту ----3.
В описанной модели не учитываются зависимости коэффициентов смкртности от плотности
населения и от других факторов, а так же не учитываются возможные и реально интересные
зависимости , которые , однако , несложно ввести в приведённую примерную схему алгоритма.
Литература
1.Гоноровский И .С.Радиотехнические цепи и сигналы .
учебник для вузов –4-е изд .-М :Радио и связь,1986.
2.Дьяконов В.П .Справочник по расчётам на микрокалькулятарах .
-3-е изд . –М : Наука,1989.
3. Дьяконов В.П .Справочник по алгоритмам и программам на языке БЕЙСИК для персональных
ЭВМ : справочник –М :наука,1987.
4.Трохименко Я . К., Любич Ф.Д .Радиотехнические расчёты на программируемых
микрокалькуляторах : справочник . –2-е изд . –М :Радио и связь, 1988.
5.Справочник по специальным функциям /под редакцией М .Абромовица и И .Стиган. перев. с англ.
–М :наука ,1979 /стр 692/.
6.Шелест А .Е . Микрокалькуляторы в физике . –М :наука,1988 /гл. 6/.
7.Инженерные расчёты на ЭВМ :справочное пособие /под ред.В.С.Троицкого –Л:Маш -е ,1979.
/параг.5.8/
8.Дьяконов В.П .Справочник по расчётам на микрокалькулятарах .
-3-е изд . –М : Наука,1989. /параг.5.11/.
9.Ильин В.Н.Основы автоматизированного схемотехнического проектирования -М :Энергия,1979.
10.Калиткин Н.Н.Численные методы. –М :наука,1978. /гл. 8/
11.Информатика и образование № 5,1990.
12.Соболь И .М .Метод Монте - Карло . –М :наука ,1972 изд . 2-е .
13.Советов Б.Я..,Яковлев С.А. Моделирование систем. Лабораторный практикум
–
М :ВШ .,1989.
JIIttN
IJ
D
,...,2,1,)(
)1(
2/
++=⋅
+−
)(
)1(
2/
:)1(
1
tN
IJ
D
N
J
It
⋅
+−
=
∑
+=
.
м
и
величин
τ
14.Стрелянов А.И .Производство вычислений на программируемых микрокалькуляторах /МК-52,
МК-54, МК-61/ -Л: Машиностроение ,1990.
15.ЭВМ в курсе общей физики /П .Е .Булкин и др. под ред.А.Н.Матвеева –М : изд -во МГУ ,1982.
ПРИЛОЖЕНИЕ 1
Разностные формулы численного решения диф .уравнений.
- -- - Устойчивость разностного метода.
Если ошибка от шага к шагу увеличивается – это приводит к неверным результатам. В этом
случае говорят, что данная задача не устойчива.Обозначим погрешность вычислений на к-м
шаге
Чтобы численный метод был устойчив необходимо , выполнения условия :
.1,)(),(
)),((),(
:1.
2
)(,
2
)(
:/,/
11
11
,1
11
2
22
2
22
порядкагоЭйлераформулаtthhOythfyy
ytxf
dt
dy
илиytf
dt
dy
порядкагоуравнениядифрешениячисленногоформуларазностнаяПростейшая
h
yyy
t
y
tyytth
h
yy
t
ty
dtyddtdyхпроизводныияапроксимацРазностая
kkkkkk
kkkk
kkkk
kkk
−−−=++=
==
−−−−
+−
=
∆
∆
=−=
−
=
∆
∆
−−−
++
−+
−
−+
... 1 переходамножителемназываютggТогда ккk
−
=
+
ε
ε
ε
Для уравнения dy/dt=f(t,y) , решаемого по формуле Эйлера первого порядка, множитель
перехода равен [6] :
Для системы диф .ур - ий :
Разностная схема Эйлера даёт :
yk+1=yk+hf1(tk,yk,zk)
zk+1=zk+hf2(tk,yk,zk)
Для того ,чтобы схема расчётов была устойчива необходимо , чтобы собственные значения
матрицы
Справка некоторых разностных формул
Уравнение 1-го порядка dy/dt=f(t,y)
1.Метод Эйлера 1-го порядка :yk+1=yk+hf(tk,yk)
2.Метод Эйера 2-го порядка (двухшаговые методы).
yk+0.5=yk+(h/2)f(tk,yk)
y
k+1
=y
k
+hf(t
k
+(h/2),y
k+0.5
)
1. 1 ≤< + g иликk εε
д y д f
h
итерацийhшагавеличинудопустимуюомаксимальнопределитьлегкооткуда
еслиустойчивМетод
h
д y
д f
h
д y
д f
g
/
2
:
111:
ˆ
,1
ˆ
1
−≤
≤+≤−
⋅=+=⋅+=
α
αα
),,(,),,( 21 zytf
dt
dz
zytf
dt
dy
==
⋅+⋅
⋅⋅+
=
⋅⋅
⋅⋅
=
−
=+=
h
д z
д f
h
д y
д f
h
д z
д f
h
д y
д f
G
Поэтому
h
дz
дf
h
дy
дf
h
дz
дf
h
дy
дf
Aматрица
единичнаяIгдеAIGпереходаматрицейсяопределяеттьУстойчивос
22
11
22
11
1
1
:
,
10
01
:
1)1)(1(42
:)1(1
ˆ
,1
ˆ
:
0...)det(
:
.]1,1[
222
1
21
011
≤−−+−+±+
⋅+=⋅+=
=++++=−
−
−
−
BAhABBABA
Прешениятиустойчивосусловиеполучимh
д z
д f
Bh
д y
д f
AОбозначив
aaaaIA
многочленатическогохарактерискорни
называютматрицызначениямимиСобственныинтервалапределызавыходилинеG
mm
mm λλλλ
λ
h
д y
д f
gh О ⋅=+== ααε ,1,)(
или : (усовершенствованный метод Эйлера-Коши) :
y*k+1=yk+hf(tk,yk)
yk+1=yk+(h/2)f((tk,yk)+f(tk+1,y*k+1))
или : y*k+1=yk+2hf(tk,yk)
yk+1=yk+(h/2)[f(tk,yk)+f(tk,y*k+1)]
3.Неявный метод Эйлера.
yk+1=yk+(h/2)[f(tk,yk)+f(tk,yk+1)] ,
Метод абсолютно устойчив, но формула не всегда разрешима в явном виде.
3.А .Общая формула сочетания явного и неявного методов.
yk+1=yk+h[А f(tk,yk)+(1-А )f(tk,yk+1)] ,
А – численный коэффициент. При А=1 получаем явный метод Эйлера. При А=0 неявный.При
0<А <1 – комбинированный [9]. При этом неустойчивость отсутствует, если 0<А <0.5 .
Примеры, когда формула разрешима в явном виде.
Пример П1.
Пример П2.Возействие сигнала произвольной формы на интегрирующую цепь .
yk+1=yk+(h/RC)[(Uk-yk)A+(1-A)(Uk-yk+1)]
4.Методы Рунге -Кутта .
Методы Рунге -Кутта имеют точность заметно выше предидущих.
Для диф .ур - ия 1-го порядка : dy/dt=f(t,y)
Формула Рунге -Кутта 2-го порядка :
д yд f
hghO
/
2
,
2
1,)(
2
2
−≤++==
α
αε
)(
2
hO= ε
.)(5.0
2
hОьпогрешностАПри =≈ ε
kk
kk
kk
U
h
h
Uоткуда
UUh
UU
U
dt
dU
+
−
=
−−+=−=
+
+
+
τ
τ
τττ
2
2
:
)(
2
1
1
1
R
U(t)
C
Uc(t
)
RC
tUctU
dt
dUc )()(
−
=
:)(
ˆ
: получимtUcyобозначив
=
)1(1
)1()([
1
A
RC
h
UAyUA
RC
h
y
y
kkkk
k
−+
−+−+
= +
Формула Рунге -Кутта 4-го порядка ,используемая в большинстве стандартных программ ЭВМ .
yk+1=yk+h/6(k1+2k2+2k3+k4)
K1=f(tk,yk) , k2=f(tk+h/2,yk+(h/2)K1 ) , k3=f(tk+h/2,yk+(h/2)K2 ) , k4=f(tk+h/2,yk+(h/2)K3 )
Для системы двух уравнений : dy/dt=f1(t,y,z) , dz/dt=f2(t,y,z) формула Рунге - Кутта 4-го порядка
имеет вид :
yk+1=yk+h/6(k1+2k2+2k3+k4)
zk+1=zk+h/6(Q1+2Q2+2Q3+Q4)
K1=f1(tk,yk,zk) Q1=f2(tk,yk,zk)
k2=f1(tk+h/2,yk+(h/2)K1,zk+(h/2)Q1 ) Q2=f2(tk+h/2,yk+(h/2)K1,zk+(h/2)Q1 )
k3=f1(tk+h/2,yk+(h/2)K2,zk+(h/2)Q2 ) Q3=f2(tk+h/2,yk+(h/2)K2,zk+(h/2)Q2 )
k4=f1(tk+h,yk+hK3,zk+hQ3 ) Q4=f2(tk+h,yk+hK3,zk+hQ3 )
Формула Рунге -Кутта для диф .ур - ия 2го порядка 5 :
Формула Рунге -Кутта для диф .ур - ия 2-го порядка 5 :
Заметим, что именно формула 4-го порядка точности , записанная дя системы произвольного
числа диф .ур - ий 1-го порядка , лежит в основе большинства стандартных программ для ЭВМ .
),(,)],(),([
2
:5.0
),(,)
2
,
2
(
:1
.5.01
),(,)]
2
,
2
(),()1[(
1
1
1
kkkkkkkkkk
kkkkkkkk
kkkkkkkkkk
ytfhfyhtfytf
h
yy
Для
ytfff
h
y
h
tfhyy
Для
либоберутобычно
ytfff
h
y
h
tfytfhyy
=++++=
=
=++⋅+=
=
==
=
⋅
+
⋅
+⋅+−+=
+
+
+
α
α
αα
αα
αα
),,( yyxfy
′
=
′
′
)()](
6
1
[
5
3211 hOkkkyhyy nnn ++++
′
+= +
)22(
6
1
43211 kkkkyy nn ++++
′
=
′
+
),,(1 nnn yyxfhk
′
⋅
=
)
2
,
8
2
,
2
(
1
12
k
yk
h
y
h
y
h
xfhk nnnn +
′
+
′
++⋅=
)
2
,
8
2
,
2
(
2
23
k
yk
h
y
h
y
h
xfhk nnnn +
′
+
′
++⋅=
),
2
,( 334 kyk
h
yhyhxfhk nnnn +
′
+
′
++⋅=
),( yxfy
=
′
′
)()]2(
6
1
[
4
211 hOkkyhyy nnn +++
′
+= +
)4(
6
1
3211 kkkyy nn +++
′
=
′
+
),(1 nn yxfhk
⋅
=
)
8
2
,
2
( 12 k
h
y
h
y
h
xfhk nnn +
′
++⋅=
)
,
(
2
3
k
h
y
h
y
h
x
f
h
k
n
n
n
+
′
+
+
⋅
=
Достоинства схем Рунге -Кутта :
1.Формулы 2-,3- и 4-го порядка точности имеют хорошую точность.
Какими из формул Рунге -Кутта целесообразно пользоваться ?
Если правая часть диф .уравнения непрерывна и ограничена вместе со своими производными
четвёртыми (и они не слишком велики ), то хорошие результаты даёт схема 4-го порядка. Если же
правая часть не имеет указанных призводных , например, для двукратно непрерывно
дифференцируемых правых частей, тогда не худшие результаты дают простые схемы меньшего
порядка точности .
2.Формулы Рунге - Кутта являются явными ,т.е . значение вычисляется по ранее найденным
значениям.
3.Все формулы допускают расчёт переменным шагом. Нетрудно уменьшить шаг там, где функция
быстро меняется и увеличить его в обратном случае . Шаг следует выбирать настолько малым,
чтобы обеспечить требуемую точность расчёта , других ограничивающих шаг условий в методе
Рунге - Кутта нет.
Так как оценка точности схем Рунге - Кутта очень громоздка, то удобнне уменьшить шаг в двое
дать апосториорную оценку точности [10,пар .8].
4.Для начала расчёта достаточно выбрать шаг h и задать t0 и y0 . Знание предыдущих значений
функции не требуется,т.е . все фо-лы Рунге - Кутта самостартующие .
Встречаются задачи , в которых функции являются достаточно гладкими , но настолько быстро
меняющимися, что схемы Рунге -Кутта , как 2,3 и 4-го порядков точности требуют неприемлимо
малого шага для получения удовлетворительного результата . Такие задачи требуют специальных
методов решения [10,пар .8].
Составитель: Будко В . Н .
Редактор: