Введение в компьютерное моделирование. Л. В. Горчаков в в в в е е д д е е н н и и е е в в к к о о м м п п ь ь ю ю т т е е р р н н о о е е м м о о д д е е л л и и р р о о в в а а н н и и е е учебное пособие
 Скачать 1.03 Mb.
|
|
Задачи для самостоятельного исследования 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 21 Рис. 3. Траектории p a cx bx x y e e k при 4, 2, 1, 3 a b c p Интегрируя по х, получаем lg lg lg a y by p x cx k или p cx a by x e y e k (14) Так как k не зависит от времени, то 0 0 0 0 p a cx by k x y e e Следовательно, мы нашли уравнения для траектории, соответ- ствующей данному начальному положению. Функция p cx x e и a by y e могут быть выражены графически (рис. 4). Рис. 4. p cx f x x e при c = 1, p = 3 22 Всякое значение за исключением экстремального она принимает 2 раза. Следовательно, если в (14) зафиксировать х, то ему соответ- ствуют два значения у и наоборот. Максимальное и минимальное значение принимаются переменной у при x p c и переменной x при y a b . Причем не существует точки перегиба. Из выражения (7) следует, что 0 dx dt , если и только если y a b . Следова- тельно, в нижней части траектории x увеличивается, а на верхней части уменьшается, таким образом, движение идет против часовой стрелки. Вычислим средние значения х и у. Поскольку движение циклическое, можно взять среднее по од- ному циклу. Пусть T – длина цикла. Из (7) имеем 1 x dx dt a by 0 0 1 T T x dx dt a by dt ; 0 lg lg 0 T x T x aT b ydt T a by (15) Но (0, T) есть полный цикл, следовательно, y a b . Аналогично x p c . Эти средние не зависят от начального положения и сов- падают с координатами точки равновесия E. Для явного получения траектории движения используем Mathcad для того, чтобы прямо решить систему уравнений (7) при следующих начальных условиях 4, 2, 1, 3 a b c p |