Файл: Л. В. Горчаков в в в в е е д д е е н н и и е е в в к к о о м м п п ь ь ю ю т т е е р р н н о о е е м м о о д д е е л л и и р р о о в в а а н н и и е е учебное пособие.pdf
ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 11.01.2024
Просмотров: 122
Скачиваний: 1
ВНИМАНИЕ! Если данный файл нарушает Ваши авторские права, то обязательно сообщите нам.
12
Рис. 1. График точного и приближенного решения
Предположение, что g(x) постоянна на отрезке от
n
x
до
1
n
x
,
означает, что скорость изменения функции у на отрезке от
1
n
x
до
n
x
постоянна и наклон касательной вычисляется в начальной точке
1
n
x
отрезка. В случае если наклон касательной меняется на неко- тором отрезке, появляется отклонение от точного решения
. Его можно уменьшить, уменьшая dx. Приближению Эйлера и истинной функции соответствует прямая и кривая рисунка. Запишем алго- ритм решения задачи для следующих начальных условий
0
max
83,
22,
0,1,
2,
0,1
s
T
T
R
t
t
В этом случае общее число шагов (точек) max
200
n
t
t
Тогда программа может иметь следующий вид
10 T=0:DT=0.1:R=0.1:N=2/DT:T1=83:TS=22 20 FOR I=1 TO N
30 CH=-R
(T1-TS)
40 T1=T1+CH
DT
50 T=T+DT
60 PRINT T,T1 70 NEXT
13
Задачи для самостоятельного исследования
1. Пусть начальная разность температур
0 61
s
T
T
. Сколько времени нужно, чтобы разность стала равной
61 2, 61 4, 61 8
?
2. Аналитическое решение уравнения (1) имеет вид
0
exp
s
s
T t
T
T
T
rt
, где
0 0
0
s
s
T t
T
T
T
T
, а
s
T t
T
Вычислите Т в момент
1
t
с шагом
0, 05; 0, 025; 0, 014; 0, 005
t
Выберите реальное r.
Постройте таблицу, содержащую разность между точными и численными решениями уравнения как функцию от
t. Будет ли эта разность убывающей функцией от
t? Если шаг уменьшить в два раза, как изменится разность? Нарисуйте график разности как функцию от
t. Если точки приблизительно расположены на убы- вающей прямой, разность пропорциональна
t (при
1
t
). Если она пропорциональна
n
t
, то численный метод называется ме- тодом n-го порядка. Какой порядок точности у метода Эйлера?
3. Какой нужно выбрать шаг
t, чтобы достигнуть точности
0,1% в момент t = 1? Какой должен быть
t, чтобы точность 0,1% была в точке t = 5?
4. Уравнение
RdQ dt
V
Q C
описывает заряд конденса- тора в RC цепи с напряжением V, t меряется в сек, R = 2000 Ом,
C = 10
–6
Ф и V = 10 В. Начальное условие Q = 0 при t = 0. Будет ли расти Q (t) со временем, увеличивается ли Q до
или будет насы- щение? Напишите программу численного решения методом Эйле- ра. Какой надо взять шаг
t, чтобы получить три правильных знака в момент t = 0,005 с. Каковы особенности численного решения при
t = 0,005, 0,004 и 0,003? Приводит ли малое изменение шага
t к большому изменению Q? Устойчив ли метод для любого шага?
5. Найдите время, необходимое для того, чтобы разность тем- ператур между T кофе и T комнаты составила
1 e
0,37 от началь- ной разности. Это время остывания называется временем релакса-
14 ции. Возьмите разные r и определите качественную зависимость времени релаксации от r.
2.1.3. Использование Mathcadа для решения задачи
об остывании чашки кофе
Оболочка Mathcad [4] позволяет легко решать численно диффе- ренциальные уравнения различного вида. При решении дифферен- циального уравнения искомой величиной является функция. Для обыкновенного дифференциального уравнения неизвестная функ- ция является функцией одной переменной. Mathcad имеет ряд встроенных функций, предназначенных для численного решения
ОДУ. В результате решения получается матрица, содержащая зна- чение функции, вычисленное на некотором множестве точек. Для поиска решения необходимо задать следующие величины:
1) начальные условия;
2) набор точек, в которых нужно найти решение;
3) само дифференциальное уравнение, записанное в определен- ном виде.
Для поиска решения можно использовать встроенную функцию
rkfixed, которая имеет следующие аргументы:
1 2
( , ,
, npoints,
)
rkfixed y x x
D
, где y – вектор начальных условий размерности n, где n – порядок дифференциального уравнения или число уравнений в системе.
Для DУ первого порядка вектор начальных условий вырождается в одну точку
0 1
y
y x
1 2
,
x x
– граничные точки интервала, на котором ищется реше- ние дифференциального уравнения. Начальные условия, заданные в векторе y, – это значение решения в точке
1
x
, npoints – число то- чек (не считая начальной точки), в которых ищется приближенное решение. При помощи этого аргумента определяется число строк
(1+ npoints) в матрице, возвращаемой функцией rkfixed. D(x, y) – функция, возвращающая значение в виде вектора из n элементов, содержащих первые производные неизвестных функций. Для полу- чения функции D (x, y) нужно разрешить уравнение относительно первой производной y
(x).
15
В результате решения получается матрица, имеющая два следу- ющих столбца:
– первый столбец содержит точки, в которых ищется решение дифференциального уравнения;
– второй столбец содержит значения найденного решения в со- ответствующих точках.
Для визуализации решения строят по этим точкам график.
Запишем уравнение (1) в виде
y
r y q
и на рабочем поле Mathcad набираем следующие данные
0
: 83 : 0,1 : 22
y
r
q
0
,
:
D x y
r
y
q
:
, 0, 2, 200,
Z
rkfixed y
D
: 0..
1
i
rows Z
Напомним, что знак (:=) набирается как (:), знак (..) набирается как (;), нижний индекс в
0
y
набирается с помощью символа ([ ).
Смысл выражений понятен из верхнего описания. Mathcad строит одну точку графика для каждого значения дискретного аргумента i, задающего график. При определении интервала изменения дис- кретного параметра i используется встроенная функция rows, кото- рая определяет количество строк в матрице Z.
Для построения графика щелкаем мышью на том месте, где его нужно разместить, и выбираем пункт ”Декартов график” из меню
“Графика”. В появляющемся пустом графике в средние поля вво- дим величины, которые нужно отобразить на графике. По оси абс- цисс – i-тую компоненту вектора
0
i
Z
, а по оси ординат – i-тую компоненту вектора
1
i
Z
. Для ввода верхнего индекса использу- ем комбинацию клавиш
16
Рис. 2. Решение c помощью Mathcad задачи о лисах и зайцах
Для рисования нескольких кривых (функций) на одном графике необходимо задать несколько выражений на оси ординат (возмож- но абсцисс). Чтобы представить графически несколько выражений на оси ординат относительно одного выражения по оси абсцисс, введите первое выражение, сопровождаемое запятой. Непосред- ственно под первым выражением появится пустое место для ввода второго выражения. Можно строить несколько графиков с парны- ми значениями функций и аргументов.
Для изменения размеров графика щелкните мышью снаружи графической области. Это закрепляет один из углов выделяющего прямоугольника. Вытяните из него прямоугольник и охватите гра- фик, затем отпустите кнопку. Цепляя за границы графика, его мож- но растянуть в нужном направлении. Отпустив кнопку, щелкните вне графика, чтобы отменить выделение.
Для изменения масштаба на графике поместите курсор в область графика, нажмите мышь, чтобы заключить график в выделяющий прямоугольник. Меню «Х–У–График» заменяет меню «Графика».
Выберите режим «Лупа» из меню – появится диалоговое окно
«Х–У–Лупа». На чертеже нажмите мышь в одном углу области, которую нужно увеличить. Удерживая кнопку, перемещаем мышь,
17 растягивая выделяющий прямоугольник. Когда вся область, кото- рую нужно увеличить, попадет в прямоугольник, отпустите кноп- ку. Координаты выбранной области отображаются в полях Мин и
Мах диалогового окна. Нажмите кнопку «Увеличить», чтобы пере- рисовать график. Чтобы сделать границы постоянными, нажмите кнопку «Принять».
Чтобы увидеть и определить координаты некоторой точки гра- фика, щелкните мышью на графике, чтобы выделить его. Выберите режим «Графики» из меню «Х–У–График», чтобы появилось окно
«Графики». Внутри чертежа нажмите кнопку мыши и переместите указатель на искомую точку. В окошках «Х-значение» и
«У-значение» отображаются значения координат этой точки. Что- бы скопировать координаты в буфер обмена, нажмите на «Копиро- вать Х» или «Копировать У». Затем можно вставлять эти значения в математическую область. Закройте окно, щелкая по кнопке си- стемного меню. Перекрестие остается на графике до тех пор, пока вы не щелкните где-либо вне графика.
2.2. Математическое моделирование в экологии
2.2.1. Моделирование нелинейных процессов в экологии
Большинство явлений природы по сути своей нелинейны [1].
Примеры нелинейных процессов – модели погоды, турбулентный режим движения жидкости. Основные понятия таких процессов легко объяснить на задаче теоретической экологии. Многие биоло- гические популяции состоят из одного поколения, которое не пере- крывается ни с предыдущим, ни с последующим. Например, энце- фалитный клещ откладывает личинки весной, а на следующую весну выводятся новые клещи. Так как процесс развития дискре- тен, то более уместно описывать его конечно-разностным уравне- нием. Простая модель развития популяции, не зависящая от плот- ности и связывающая численность популяции в
1
n
– поколении с предыдущим n-поколением, записывается в виде
1
P
P
n
n
a
, где
а – константа. При
1
a
получается геометрический рост популя- ции, что нереально. В более реалистической модели численность популяции ограничивается окружающей средой. Простая дискрет-
18 ная модель прироста популяции, зависящая от плотности, записы- вается в виде
1
P
P
n
n
n
a bP
. Второй член представляет уменьшение естественного прироста за счет перенаселения или распространения болезней. Положим
n
n
P
a b x
и получим уравнение
1 1
n
n
n
y
ax
x
. Введем параметр роста
4
r
a
и получим уравнение вида
1
n
n
x
f x
, где
4 1
f x
rx
x
. Ес- ли
1
n
x
, то значение
1
n
x
будет < 0. Чтобы избежать этой нере- альной ситуации, наложим условия ограничения изменения пере- менной x и параметра r отрезками
0 1
x
и
0 1
r
Постройте решение уравнения для начальных условий
0 0, 75
x
и значений
0,1;0,8;0,9
r
в зависимости от n с помо- щью программы и с помощью Mathcad.
2.2.2. Борьба популяций в теоретической экологии
Рассмотрим взаимодействие двух видов животных, один из ко- торых служит пищей для другого [5]. Пусть животные первого ви- да кролики, второго – лисы. Аналогичной задачей является распро- странение эпидемий с помощью бактерий, здесь можно рассматри- вать население как животные первого вида, а бактерии второго вида.
Рассмотрим животных первого вида. Если нет лис, то популяция кроликов будет увеличиваться (другие факторы, влияющие на рост этих видов животных, остаются неизменными, как переменная ве- личина учитывается только число животных). Простейшим пред- положением является, что скорость роста популяции кроликов пропорциональна размеру популяции. Пусть х – число кроликов, у
– число лис в момент времени t. Тогда, если
0
y
, то
,
0
dx dt
ax a
(5)
С другой стороны, не имея пищи, лисы будут вымирать, так что получаем, если
0
x
, то
,
0
dy dt
py p
(6)
Если же имеются и кролики и лисы, то необходимо принять во внимание их взаимодействие. Предположим, что число съеденных кроликов пропорционально величине xy. Тогда необходимо доба-
19 вить в (6) член, пропорциональный ху, чтобы позволить увеличение числа лис при наличии пищи и вычесть такой же член из (5), чтобы учесть съедаемых кроликов
,
,
0
dx dt
ax bxy
a b
,
,
0
dy dt
cxy
py
c p
(7)
Мы получили динамическую модель взаимодействия двух видов животных, она описывается системой дифференциальных уравне- ний первого порядка.
Для решения такой системы могут применяться различные ме- тоды. Общий вид полученной системы уравнений можно записать так
,
dx dt
F x y
,
dy dt
G x y
(8)
Подсчитаем количество животных каждого вида в данный мо- мент времени, который можно принять за
0
t
, и получим
0 0
x
x
кроликов и
0 0
y
y
лис. Нашей главной задачей яв- ляется определение численности обеих популяций в будущем. При делении уравнений друг на друга время исключается и получаем уравнение вида
,
,
dy dx
G x y
F x y
(9)
Для таких систем уравнений доказаны две теоремы.
Теорема 1. Если в окрестности точки плоскости
0 0
,
x y
част- ные производные функций F и G непрерывны, то существует един- ственное решение, проходящее через
0 0
,
x y
при
0
t
. Решения либо не зависят от времени, либо описываются гладкой кривой.
Кроме того, решение
x t
и
y t
непрерывно зависят от началь- ного положения. Существуют три типа решений: устойчивые, не- устойчивые и циклические.
Теорема 2. Поведение траектории вблизи точки равновесия можно определить, рассматривая только линейные члены разложе- ния функций F и G в ряды Тейлора в точке равновесия (которая случается при
0
F
G
). Решения получаемых линейных урав-
20 нений в окрестности точки равновесия имеют то же поведение, что и точные решения.
Все сказанное справедливо и для n уравнений с n неизвестными.
1 2
,
,...,
,
1..
i
i
n
dx dt
F x x
x
i
n
(10)
Интересной точкой равновесия является точка
,
E
p c a b
. Ис- следуем поведение вблизи E. Положим
u
x
p c
и v
y a b
Тогда
v v
v
du dt
dx dt
u
p c b
d
dt
dy dt
cu
a b
и
Линейные части этих уравнений есть
v, v
du dt
bp c
d dt
ac b u
(11)
Рассмотрим эти уравнения как точные (согласно теореме 2). Диф- ференцируем первое уравнение и, подставляя v
d dt
из второго, получим
2 2
d u dt
apu
(12)
Следовательно, движение будет периодическим. Поскольку значе- ние начального момента неважно, начнем с момента, когда
0
u
Решение (12) тогда есть
sin
u
A
apt
и из (7) получим
v cos
B
apt
. Следовательно,
2 2
2 2
v
1
u
A
B
, т.е. траекто- рия есть эллипс (рис. 3).
Таким образом, вблизи точки равновесия Е траектории суть пе- риодические движения вокруг точки равновесия. В первом при- ближении эти траектории являются эллипсами с периодом обраще- ния
2
ap
Для нахождения точных траекторий образуем уравнение со- гласно (9)
dy dx
y cx
p
x a by
(13)
Отсюда
0
a by
y dy dx
p cx
x