Этюд
5
Решение дифференциальных
уравнений
Рис. 5.1. Задача об эпидемии –
разностная схема
Рис. 5.2. Задача об эпидемии –
решение системы дифференциальных уравнений
Рис. 5.3. Решение краевой задачи об
эпидемии разностной схемой
Рис. 5.4. Решение краевой задачи об
эпидемии функцией rkfixed
Рис. 5.5. Решение краевой задачи об
эпидемии функцией sbval
Рис. 5.6. Моделирование развития
финансовой пирамиды
Рис. 5.7. Развитие финансовой
пирамиды
Одно из основных потребительских качеств
компьютера – это отношение цены к производительности. С подобным критерием
можно подойти и к примерам, входящим в пакеты программ. Только вместо
производительности нужно рассматривать занимательность задачи, а вместо цены –
размер соответствующего файла. В этом смысле рекордсмен пакета Mathcad– задача «Эпидемия».
Ее уникальность еще и в том, что она – единственная, автор которой (John Truxal)
отмечен в разделе Acknowledgments документации пакета Mathcad.
Суть задачи. В городе, население которого
составляет 20 000 человек (см. пункт 1 на рис. 5.1), появляются 50 инфекционных
больных, что вызывает эпидемию. Предположим, что прирост больных за день
пропорционален произведению числа здоровых (еще не переболевших и не
приобретших иммунитет) на число больных. Коэффициент пропорциональности Пр интегрирует разного рода меры профилактики. Если, к примеру, жители
города будут носить марлевые повязки или сделают прививки, то этот коэффициент
уменьшится. Спрашивается, как развивается эпидемия: как изо дня в день меняется
число больных.
Все это словесное описание модели эпидемии
легко вмещается в Mathcad-документ (рис. 5.1) с двумя формулами (пункт 2),
объединенными в вектор, что эквивалентно BASIC-конструкции:
For t = 1 To 13
Больные(t + 1) = Пр * Больные(t) * Здоровые(t)
Здоровые(t + 1) = Здоровые(t) - Больные(t)
Next
Если бы два выражения в пункте 2 на рис. 5.1
не были заключены в скобки, то их выполнение сразу бы прерывалось сообщением об
ошибке. Система Mathcad пыталась бы сначала полностью заполнить вектор Больные, а уже потом – вектор Здоровые. Скобки изменяют порядок
счета: он ведется не по строкам, а по столбцам: сначала заполняются вторые
элементы векторов Больные и Здоровые (первые элементы заполняются в пункте 1 – начальные условия), а потом
третьи и т.д. Этими скобками мы меняем естественный порядок выполнения
операторов[1]
¾ они выполняются не слева
направо и не сверху вниз, а крест накрест.
Результаты расчета графически отображены[2]
в пункте 3. На тринадцатый день (спад эпидемии) в городе было 105 больных.
Критическая точка – девятый день (3972 больных), ради поиска которой и затевают
весь этот расчетный сыр-бор: моделируя эпидемию, мы можем распланировать работу
санитарных служб города – подвезти в аптеки лекарства, отозвать врачей из
отпуска, выписать из больниц выздоравливающих и т.д. В оригинальном решении
задачи дополнительно прослеживается динамика изменения числа умерших во время
эпидемии. Но мы исключим эту печальную кривую. Тем более что она и
математически не вписывается в задачу: динамика числа умерших просчитывается
отдельно от динамики больных и здоровых.
Описанная задача об эпидемии сводится к
решению задачи Коши для системы двух обыкновенных
дифференциальных уравнений (ОДУ) первого порядка относительно двух
неизвестных функций Больные(t) и Здоровые(t).
Больные’(t) = Больные(t) × (Пр × Здоровые(t)-1)
Здоровые(t) = -Пр × Больные(t) × Здоровые(t)
В пункте 2 на рис. 5.1, по сути, реализован метод Эйлера[3] с
единичным шагом интегрирования Dt.
Больные(t+1)
= Больныеt + Dt × Больные¢t
Здоровые(t+1)
= Здоровыеt + Dt × Здоровые¢t
Заглядывая вперед, в этюд 6 (вернее, готовя
читателя к нему), покажем, как вектор-выражение пункта 2 на рис. 5.1
реализуется средствами программирования (цикл с параметром):
for t Î 1..13
|
Больныеt+1
¬
Пр×Здоровыеt
×Больныеt |
|
3доровыеt+1
¬
3доровыеt - Пр×Здоровыеt ×Больныеt |
В среде Mathcad до версий PLUS 5.0
дифференциальные уравнения без особых ухищрений можно было решать только методом
Эйлера, который имеет низкую точность и производительность (плата за простоту).
Инструментарий для решения дифференциальных уравнений (систем) различного
порядка и различными методами появился в арсенале Mathcad PLUS 6.0. В него
входят 13 встроенных функций (Bustoer,
bustoer, bvalfit, multigird, relax, Rkadapt, rkadapt, rkfixed, sbval, Stiffb,
stiffb, Stiffr
и stiffr), работа одной из которых (rkfixed – метод Рунге ¾ Кутты (rk) четвертого порядка с фиксированным (fixed) шагом интегрирования) показана на рис. 5.2. У этой функции
пять аргументов: вектор начальных значений искомых решений, абсцисса начальной
точки интегрирования, абсцисса конечной точки интегрирования, число шагов
интегрирования и функция-вектор правых частей системы. Функция rkfixed возвращает в матрицу (у нас она имеет имя Z) с P+1 столбцами и n строками (P – количество уравнений или порядок уравнения) таблицу решений системы:
первый (вернее, нулевой) столбец – это значения аргумента t (их задает пользователь через величины tнач, tкон[4] и n), а последующие столбцы – значения ординат решения.
В пункте 5 рис. 5.2 столбцы матрицы Z (время, число больных и число здоровых) разнесены по отдельным векторам
и отображены графически, что позволяет проследить динамику развития эпидемии.
На рис. 5.2 получены несколько иные
результаты, чем на рис. 5.1, хотя характер кривых сохранился: максимум больных
(3119) наблюдается не на 9-й, а на 7-й день, на 13-й день мы имеем не 105, а
209 больных. Это объясняется и различными значениями точности расчетов (на рис.
5.1 делалось 13 шагов интегрирования, а на рис 5.2 ¾ 500) и различными
примененными методиками (Эйлер против Рунге и Кутта[5]).
В функцию rkfixed заложен широко
распространенный метод решения дифференциальных уравнений – метод Рунге ¾ Кутта[6]. Несмотря на то что это не
самый быстрый метод, функция rkfixed почти всегда справляется с
поставленной задачей. Однако есть случаи, когда лучше использовать более сложные
методы. Эти случаи попадают под три широкие категории: система может быть жесткой[7]
(Stiffb, Stiffr), функции системы могут быть гладкими
(Bulstoer) или плавными (Rkadapt). Нередко приходится пробовать
на одном дифференциальном уравнении (одной системе) несколько методов, чтобы
определить, какой метод лучше (быстрее, точнее). Примерно так мы сравнивали в
этюде 3 разные способы поиска оптимумов функции.
Когда известно, что решение гладкое, используется
функция Bulstoer, куда заложен метод Булирша ¾ Штёра, а не Рунге ¾ Кутты, используемый
функцией rkfixed. В этом случае решение
будет точнее. Список аргументов и матрица, получаемая при работе с функцией Bulstoer, такие же, как и при работе
с rkfixed.
Можно решить задачу более точно (более
быстро), если уменьшать шаг (у нас это Dt) там, где производная
меняется быстро, и увеличивать шаг там, где она ведет себя более спокойно. Для
этого предусмотрена функция Rkadapt (adaption –адаптация). Но, несмотря на то что при решении
дифференциального уравнения функция Rkadapt использует непостоянный
шаг, она тем не менее представит ответ для точек, находящихся на одинаковом
расстоянии, заданном пользователем. Аргументы и матрица, возвращаемая функцией Rkadapt, такие же, как при rkfixed.
При решении жестких систем следует
использовать одну из двух встроенных функций, разработанных специально для
таких случаев: Stiffb и Stiffr.
Они используют метод Булирша ¾ Штёра (b) или Розенброка (r).
Форма матрицы-решения, полученной с помощью этих функций, идентична матрице,
полученной через rkfixed. Однако
Stiffb и Stiffr требуют
дополнительного аргумента J:
Stiffb(x,
tнач, tкон, n, f, J)
Stiffr(x,
tнач, tкон, n, f, J),
где
x – вектор n начальных значений;
tнач, tкон – конечные точки интервала, в котором должно быть найдено решение
дифференциальных уравнений. Начальные значения x определяются в точке tнач;
n – количество точек за начальной точкой, в которых
должно быть определено решение. Это определяет число рядов (1 + n) матрицы, которую генерируют функции Stiffb и Stiffr;
f(t, x) – функция-вектор правых частей системы;
J(t, x) – матрица-функция размерности n×(n+1), в которой содержится матрица Якоби правых частей дифференциальных
уравнений.
В функциях для решения дифференциальных
уравнений, описанных ранее, предполагалось, что необходимо получить таблицу
значений x(t) относительно значений t, расположенных на одинаковом расстоянии друг от друга на отрезке
интегрирования, ограниченном точками tнач и tкон. Но часто требуется найти решение только в конечной точке tкон. Хотя описанные функции определяют значение x(tкон), они проделывают при этом много ненужной работы, высчитывая
промежуточные значения x(t).
Если необходимо узнать только значение x(tкон), используются следующие функции:
bulstoer(
x, tнач, tкон, acc, f, kmax, save)
rkadapt(
x, tнач, tкон, acc, f, kmax, save)
stiffb(
x, tнач, tкон, acc, f, J, kmax, save)
stiffr
( x, tнач, tкон, acc, f, J, kmax, save),
где
x – см. выше;
tнач, tкон – конечные точки интервала, на
котором должно быть найдено решение дифференциальных уравнений. Начальные
значения t определяются в точке tнач;
acc – контролирует точность решения. При небольших
значениях acc делаются более мелкие шаги
вдоль траектории, что увеличивает точность решения;
f(t, x) – см. выше;
J(t, x) – см. выше;
kmax – максимальное количество промежуточных точек, в
которых будет найдено решение. Значение kmax устанавливает верхнюю
границу количества рядов матрицы, получаемой этими функциями;
save – наименьшее допустимое расстояние между
величинами, в которых должно быть найдено решение. Значение save устанавливает нижнюю границу разницы между любыми двумя числами в
первой колонке матрицы решения.
Вернемся к нашему примеру с эпидемией. Если
известно начальное число больных (N=50 – см. пункт 1 на рис. 5.1),
то это значит, что их сразу пересчитали. А если это так, то их знают поименно.
Но в этом случае больные должны быть изолированы, и никакой эпидемии не будет
вообще (Пр=0). В реальной задаче мы можем
знать число жителей в городе (число здоровых на день начала эпидемии) и число
больных в какой-то последующий день, когда становится ясно, что разразилась
эпидемия. Другими словами, нам что-то известно о параметрах на краях отрезка, охватывающего некий динамический
процесс.
На рис. 5.3. реализован метод
последовательных приближений для решения по разностной схеме такой краевой задачи: в начале эпидемии в
городе 20 000 здоровых жителей, а в конце эпидемии – 100 больных. Спрашивается,
сколько больных было в начале эпидемии.
Ответ (51 больной) получен за семь приемов:
задается начальное число больных, которое корректируется в зависимости от того,
какое число больных оказывается в конце эпидемии. На рис. 5.3 можно, конечно,
не дублировать Mathcad-оператры, а просто вручную подправлять первое
приближение.
С помощью функций Bustoer, bustoer, Rkadapt, rkadapt, rkfixed, Stiffb, stiffb, Stiffr и stiffr, решающих задачу Коши, краевую задачу можно также решить последовательными приближениями (см.
рис. 5.4 с функцией rkfixed):
Нижняя кривая на рис. 5.2 похожа на
траекторию полета снаряда[8],
поэтому метод последовательных приближений, приложенный к краевой задаче,
называют также методом стрельбы:
можно менять искомые начальные условия, последовательно приближаясь к решению,
имея на другом конце отрезка «корректировщика огня», в лексиконе которого три
слова: «перелет», «недолет» и «попал» («почти попал», учитывая заданную
точность расчета – см. комментарии на рис. 5.3). По правде говоря, метод
стрельбы берет свое название не от вида кривой, а от особенностей решения
краевой задачи применительно к дифференциальному уравнению второго порядка,
когда приходится менять угол наклона ствола пушки – значение первой производной
на конце отрезка интегрирования. В нашей задаче об эпидемии (система двух
обыкновенных дифференциальных уравнений) для попадания в цель приходится менять
не угол наклона пушки, а ее подъем над землей – число больных в начале
эпидемии. Тут, как правило, используют метод половинного деления, когда отрезок
«перелет-недолет» делят пополам.
Метод стрельбы заложен и во встроенную функцию
sbval, работа которой показана на
рис. 5.5. Она требует начального приближения и соблюдения некоторых
условностей, связанных с нашими знаниями состояний решаемой системы на концах
отрезка интегрирования. Функция sbval не совсем обычная. Мы
привыкли к тому, что, например, операторы A:=sin(X) и A:=sin(30) идентичны, если переменная X имеет значение 30, несмотря на то, что в первом случае у функции sin в качестве аргумента выступает переменная, а во втором – константа.
Функция же sbval в этом отношении аномальна.
Она возвращает значение в зависимости не только от конкретных значений
аргументов, но и от того, введены они в виде переменных (x0) или в виде констант (20 000). Из-за этого даже специалистам по дифференциальным уравнениям часто
приходится ломать голову, чтобы сообразить, как можно условия краевой задачи
«запихнуть» в функцию sbval. Здесь фирма MathSoft
несколько отошла от первоначальной идеологии пакета Mathcad, подразумевавшей,
что запись условия задачи на экране дисплея должна выглядеть естественно.
Типичный пример – численное решение в среде
Mathcad системы алгебраических уравнений: сначала задается начальное
приближение к корню (ведь корней может быть много), потом за ключевым словом Given (дано) в естественном виде пишется сама система уравнений, после
которой помещается функция Find, возвращающая значения
своих аргументов, превращающих ранее записанные уравнения (до слова Given) в тождества. Функция Find также аномальна: она
возвращает значения, зависящие не только от значений аргументов, но и от того,
что вблизи нее находится. Это, конечно, грубейшее нарушение правил построения
функциональных зависимостей. Но разработчикам Mathcad пришлось пойти на это
ради того, чтобы алгебраические уравнения в Mathcad-документе хранились в
естественной форме. При создании математических пакетов типа Mathcad приходится
решать не только чисто математические проблемы (разработка алгоритмов решения и
реализация их средствами вычислительной математики), но и проблемы единства формы представления задачи (интерфейс
пользователя) и методов ее решения.
Функция sbval не решает краевую задачу, а
только находит недостающие значения на краю отрезка. После этого краевая задача
переходит в задачу с начальными условиями (задача Коши), которая и решается в
пункте 7 рис. 5.5. Заодно проверяется точность решения: задали 100 больных –
получим почти 100 (ошибка во втором знаке после запятой).
Для решения краевой задачи в среде Mathcad
есть еще одна функция – bvalfit. Она используется в тех
случаях, когда нет всей необходимой информации для функции sbval, но известно решение задачи в промежуточной точке. Функция bvalfit решает задачу с двумя граничными точками, начиная с конечных точек и
следуя траекториям решения и его производным в промежуточных точках:
bvalfit( x1, x2, tнач, tкон, tf, f, load1, load2, score),
где
x1- вектор предполагаемых начальных значений, не
определенных в точке х1;
x2 – то же для точки х2;
tнач, tкон – конечные точки интервала, в котором должно быть вычислено решение
дифференциального уравнения;
tf – точка между tнач и tкон, в которой траектории решений, начатые с tнач, и траектории, начатые с tкон, должны совпадать;
f(t, x) – векторная функция, состоящая из n элементов, содержащая первые производные неизвестной функции;
load1(tнач, x1) – векторная функция, n элементов которой соответствуют значениям n неизвестных функций в точке х1. Некоторые из этих значений
будут постоянными, заданными начальными условиями. Другие будут неизвестны
вначале, но будут найдены. Если значение неизвестно, то следует использовать
соответствующее предполагаемое значение из вектора x1;
load2(tкон, x2) – аналог load1, но для значений n неизвестных функций в точке
х2;
score(tf, x) – векторная функция, состоящая
из n элементов. Эта функция применяется для того,
чтобы задать, как решения должны совпадать в точке tf. Обычно нужно определить score(tf,
x) := x,
чтобы все решения неизвестных функций совпадали в точке tf.
Функция bvalfit по своей сути решает две
краевые задачи на двух смежных отрезках интегрирования. Она особенно полезна,
когда производная имеет разрыв где-то внутри интервала интегрирования. При
развитии эпидемии (наша задача) такое может произойти, если, например, комиссия
из центра снимет с работы медицинское начальство города и тем самым изменит значение
коэффициента Пр.
Второй тип задач с граничными условиями
появляется при решении дифференциального уравнения с частными производными. В
этом случае приходится фиксировать значение решения не в двух точках, как было
сделано в предыдущей «плоской» задаче, а в совокупности точек, составляющих
границу: в углах прямоугольника, например. В среде Mathcad имеются две функции
(relax и multigird) для решения таких уравнений. Но разговор о них выходит за рамки этой
книги. Скажем лишь то, что эти функции предназначены для решения
дифференциальных уравнений конкретного вида – уравнений Лагранжа и Пуассона.
Задачи другого вида требуют составления индивидуальной схемы решения с
привлечением той же «плоской» функции rkfixed.
А теперь наши читатели получат домашнее
задание.
Даны математическая модель развития
финансовой пирамиды и ее реализация по разностной схеме. Необходимо составить и
решить систему дифференциальных уравнений, описывающих развитие эпидемии,
пардон, финансовой пирамиды.
К задачам,
решаемым на персональном компьютере, не в последнюю очередь относятся и финансовые, связанные с ведением бюджета
семьи, фирмы, завода и т.д. Итак, финансы. Компьютер может помочь учитывать,
сохранять и приумножать их минимум тремя способами:
1.
На компьютере заводится электронная версия
бухгалтерской книги, куда записываются все доходы и расходы. Для этих целей
годятся такие программы, как Excel, Lotus 1-2-3, QuattroPro
(электронные таблицы), 1C:Бухгалтерия
(специализированные пакеты[9])
и т.д. Пакет Mathcad, кстати, ведет свою родословную от электронных таблиц,
которые также можно считать «волшебной» счетной доской. Первое его отличие
состоит в том, что в среде электронных таблиц знак «равно» ставится перед
выражением, а в среде Mathcad – после. Второе отличие в том, что Mathcad
поддерживает графический режим интерфейса: =
5×25
+ 3Ù12 (Excel); 5×25 + 312 = (Mathcad). Третье отличие, максимально приближающее
стиль Mathcad к стилю написания формул на обычном листе бумаги, заключается в
том, что формулы и константы электронных таблиц жестко привязаны к сетке
столбцов и строк. Об этом мы подробно поговорим в приложении 7, когда будем
описывать методику совместной работы Mathcad и Excel.
2.
Компьютер подсоединяется к банковской
электронной сети так, что домашняя хозяйка или бухгалтер фирмы могут сами
проводить необходимые платежи, не выходя из дома. Здесь на первый план
выступает проблема безопасности:
платежное поручение должно быть зашифровано и сопровождаться электронной
подписью. Об этом мы еще поговорим в этюде 7.
3.
На компьютере можно смоделировать и
просчитать ту или иную финансовую операцию (покупка облигаций или акций,
открытие счета или взятие кредита в банке и т. д.) и посмотреть, во что это
выльется.
Среди примеров, входящих в пакет электронных
таблиц Excel, есть задача, связанная с покупкой ценных бумаг. Рассчитывается,
сколько и каких акций нужно купить, имея в запасе ограниченное количество
свободных денег, чтобы сумма будущих дивидендов была максимальной (задача
линейного программирования – см. рис. 2.9 и 2.10 в этюде 2, 3.10 и 3.11 в этюде
3).
Среди примеров, входящих в пакет
математической программы Mathcad, есть задача, связанная с моделированием
развития эпидемии. Задается начальное число здоровых и больных, а далее
просчитывается, как развивается эпидемия – как меняется число больных по дням.
Мы только что эту задачу разобрали.
Попробуем скрестить эти две задачи и получить
новую, связанную с моделированием развития финансовой
пирамиды, которая во многом похожа на эпидемию. Сделаем мы это в среде
пакета Mathcad, который вполне подходит для таких целей. Более того, к пакету
Mathcad можно прикупить электронный учебник «Personal Finance», который поможет научиться
вести финансы и из которого при необходимости в рабочий документ переносятся
нужные формулы и константы.
Автор имел опыт участия в финансовых
пирамидах. Первый опыт ¾ очень давний, очень
невинный и без особых последствий. В школьные годы кем-то в классе было
предложено не тратить выдаваемые родителями ежедневные обеденные 20 копеек, а складывать
их и отдавать все разом поочередно каждому участнику этой финансовой кампании.
Так воплощалась старая, как мир, детская мечта: «Вот если бы каждый взрослый дал мне по копеечке... Он бы этого и не
заметил, а у меня бы оказался целый капитал». Но в классе вышло так, что
участники этой складчины, получившие деньги, из игры выбывали, она потихоньку
глохла и кто-то (а среди них оказался и автор) остался с носом.
Второй опыт участия в финансовых пирамидах
автор приобрел вместе со всей страной 17 августа 1998 года.
Сейчас, слыша о крахе «Тибетов», «Светлан» и
прочих «Чар», автор вспоминает детскую мечту, школьные двугривенные и не только
их.
Банковская система любой страны, как на трех китах, покоится на трех числах. Первое число N1 – плата за кредит. Взял в банке сто рублей – будь любезен в конце года
верни 100 + N1 рублей. Второе число N2 – процент по вкладу. Положил в банк сто рублей – получи в конце года 100 + N2 рублей. Разница между первым и вторым числами (N1 > N2) заставляет банки прибыльно работать. Третье число N3, подпирающее снизу два первых и заставляющее людей нести деньги в
банк, – это величина инфляции[10].
В нормальной экономической ситуации низкий уровень инфляции и не очень высокая
плата за кредит держат в узких рамках процент по вкладу:
N1 > N2 > N3.
Если же инфляция составляет 1000 и более
процентов в год, то многие люди, забывая о ненормальности такой ситуации, легко
верят в 1000 и более процентов годовых по вкладу (ведь величина N2 должна быть больше величины N3) и «ложатся в основание»
очередной финансовой пирамиды. Если, конечно, законодательством страны
позволительно такие пирамиды строить. Есть и менее наивные люди, понимающие,
что пирамида – это особый род игры, где нужно уметь «вовремя смыться».
Итак, строим финансовую пирамиду.
В пункте 1 рис.
5.6 вышеприведенного протокола работы в среде Mathcad определяются константы –
ее имя, знак присвоения (:=) и величина. Комментарии расшифровывают их.
В пункте 2
определяется состояние пирамиды на первый день: вводятся индексные переменные –
первые значения векторов M, NK и MMM.
В пункте 3
записана динамика изменения курсов продажи и покупки акций – функции P(t) и K(t): объявляется о выпуске
акций (билетов) номиналом в 100 рублей со следующим курсом продажи P и покупки K:
Таблица 5.1
|
Число
дней, прошедших с начала эмиссии акций (билетов) |
1 |
2 |
3 |
... |
51 |
... |
365 |
... |
|
Продажа
(руб.) |
105 |
107 |
109 |
... |
205 |
... |
833 |
... |
|
Покупка
(руб.) |
100 |
102 |
104 |
... |
200 |
... |
828 |
... |
Из таблицы 5.1 видно, что купленная акция
может дать дивиденд в 723% годовых при номинальной своей цене в 100 рублей.
Если уровень инфляции достаточно высок, то люди верят в реальность таких
огромных дивидендов и пирамида растет. Но опасность краха этой затеи ощущают
почти все и отдают свои деньги не на год, а, допустим, на 50 дней (переменная Время – среднее время между покупкой и продажей акций – см. пункт 1). За
этот период по каждой акции можно «наварить» магические 100 рублей,
фигурирующие во многих пословицах и поговорках.
Далее векторы NK и NP заполняются по простой
разностной схеме: известно предыдущее значение элемента вектора (на день t) – рассчитывается его очередное значение (на день t+1).
В городе, где
строится пирамида, миллион жителей (N – см. пункт 1), среди
которых царит некий ажиотаж, подогреваемый вышеприведенной таблицей курсов.
Языком математики его можно описать формулой, связывающей число проданных
населению акций в конкретный день (NK) с общим числом проданных
акций (сумма NK за предыдущие дни) и
условным числом жителей, не купивших пока акции (N минус сумма NK за предыдущие дни). Повторяем, развитие финансовой пирамиды во многом
напоминает развитие эпидемии, когда число заболевших (купивших акции) в
конкретный день пропорционально числу больных в городе (числу проданных акций),
перемноженному на число еще не переболевших (не купивших акции). В случае
эпидемии коэффициент пропорциональности зависит от мер профилактики. В случае финансовой
пирамиды этот коэффициент (мы его условно назовем коэффициентом ажиотажа – KA) зависит от уровня
инфляции, рекламы (вспомним Марину Сергеевну, Леню Голубкова и прочих
«бабочек»), наличия других параллельных пирамид, от срока, прошедшего с момента
шумного краха предыдущей пирамиды, и т.д. Многие экономические явления
(кризисы, банкротства) прокатываются волнами. Период пика волн финансовых
пирамид составляет, по различным оценкам, от 25 до 30 лет, что связано,
во-первых, с приходом к активной жизни свежих, незатронутых пирамидами сил, и,
во-вторых, с короткой людской памятью. На таких волнах многих ждет финансовое
кораблекрушение. Другие же (а их намного меньше – и в этом фокус пирамид),
подобно отважному и ловкому серфингисту, получают «финансовое» удовлетворение.
За волной
купивших акции «катит» волна желающих их продать – вернуть свои «кровные» и
причитающиеся дивиденды. Здесь мы также до предела упростим модель и будем
считать, что волна продающих акции отстает от волны их купивших на число дней,
хранящихся в переменной Время:
NPt+1 = 0, если t £ Время
NPt+1 = NKt-Время, если t > Время.
Волны
покупателей и продавцов акций могут иметь разные формы – подчиняться, например,
нормальному закону распределения (см. в этюде 6 рис. 6. 41 в этюде 6). Главное здесь
раздвоенность волн: человек сначала
покупает акцию (билет) и только потом ее продает.
Ну а теперь
можно подсчитывать барыши и кататься на волнах финансовой пирамиды.
Несложно
вычислить, сколько денег (вектор M ¾ см. пункт 5 на рис. 5.7)
будет на счету организаторов пирамиды завтра (t
+ 1), если
известно, сколько их в наличии сегодня (t), и если известен курс
акций и количество покупок и продаж:
Mt+1 = Mt + NKt × K(t) -
NPt × P(t)
Люди, покупающие
акции, приносят деньги в кассу. Люди, акции сдающие, забирают деньги из кассы.
Но есть еще один человек, залезающий в кассу. Это – организатор пирамиды,
имеющий свой «профит», что выражается в том, что из кассы ежедневно изымаются
три процента (Доход := 0.03 – см. пункт 1) наличных
денег:
Доход × Mt
Естественно,
доход изымается, если (if) в кассе есть деньги. В
реальной жизни, конечно, касса худеет на значительно большие суммы – налоги,
оплата текущих расходов, реклама и т.д. Расход
:= 300 000 –
см. пункт 1.
В 1202 году Леонардо
Пизанский (1180-1240) описал одну из первых моделей развития замкнутой
биологической системы, населенной условными кроликами. Если соответствующим
образом определить их плодовитость и долголетие, то численность популяции
кроликов будет меняться из поколения в поколение по строгому закону:
Таблица 5.2
|
Поколение |
1 |
2 |
3 |
4 |
5 |
6 |
7 |
... |
27 |
... |
|
Число кроликов |
1 |
1 |
2 |
3 |
5 |
8 |
13 |
... |
196 418 |
... |
Читатель,
конечно, уже догадался, что речь идет о числах
Фибоначчи: Леонардо Пизанский более известен под именем Фибоначчи
(Fibonacci – сокращение от filius Bonacci – сын Боначчи). В новом поколении
кроликов их число будет равно сумме числа кроликов в двух предыдущих поколениях
(см. программу на рис. 6.11 в этюде 6). Со временем про этих условных кроликов
забыли, но числа Фибоначчи (1, 1, 2, 5, 8, 13 и т.д.) нашли применение в
прикладной математике (см. этюд 6).
Наша модель
развития пирамиды также позволяет сгенерировать некий числовой ряд,
отображающий состояние дохода организаторов этой финансовой операции (вектор
МММ):
Таблица 5.3
|
День |
Доход (округлено до
рублей) |
Примечание |
|
1 |
70 000 000 |
Начало пирамиды |
|
2 |
72 100 000 |
|
|
3 |
74 128 022 |
|
|
4 |
76 086 206 |
|
|
... |
... |
|
|
168 |
303 058 831 |
|
|
169 |
303 485 168 |
|
|
170 |
303 635 916 |
День икс |
|
171 |
303 635 916 |
|
|
272 |
303 635 916 |
|
|
... |
... |
|
Назовем числа второго столбца таблицы 5.3 числами Мавроди[11].
Будем надеяться, что со временем о финансовых пирамидах забудут, но числа
Мавроди войдут в историю. Тем более, что сам Сергей Мавроди по образованию
математик.
В пунктах 5-6
графически отображено развитие пирамиды. Просматривая матрицу M, можно определить «день
икс», когда прибыль организатора достигает максимума (у нас это 170-й день –
см. пункт 7), и когда пирамиду пора разваливать – уходить на «дно»,
баллотироваться в депутаты или уезжать за границу. Благо денег на это
«наварено» достаточно – почти триста миллионов при всего лишь трехпроцентной
норме прибыли.
Мы же никуда пока не уезжаем, остаемся у
своего компьютера и, собираясь вкладывать деньги в какое-то надежное или
сомнительное предприятие, сначала должны просчитать, что из этого может выйти.
Так мы легко можем вернуть и приумножить деньги, потраченные на приобретение
компьютера и программы Mathcad[12],
а также на операционные системы Windows, под управлением которой Mathcad
работает.
Кроме того, читатель должен не забывать о
задании, которое он получил, – решить задачу о финансовой пирамиде не по
разностной схеме, а через дифференциальные уравнения.
[1] Об этом атрибуте
программирования будет подробнее сказано в этюде 6.
[2] Здесь используются два типа
декартова графика: линия с точками в виде квадратиков и так называемая bar-diagram
– плоская столбчатая диаграмма.
[3] Или Ойлера – что в него
заложено, можно увидеть в этюде 6 на рис. 6.3.
[4] Здесь лучше написать Тнач и Ткон, опустив индекс.
[5] Многие математики полагают,
что есть только один метод Эйлера. Все остальное ¾ модификации этого метода.
[6] Или Рунге – Кутты, что в него заложено – см. рис. 6.3 в этюде 6.
[7] Термин «жесткий» происходит из механики, где численное решение некоторых систем дифференциальных уравнений требует разного шага интегрирования по разным искомым функциям.
[8] Она может быть совсем не похожа на траекторию полета снаряда, но… читаем дальше.
[9] В этюде 6 мы рассмотрим
методику составления одной функции для такого рода расчетов, возвращающую сумму
налога.
[10] Некоторые банки
Великобритании предлагают своим клиентам такие условия: процент по вкладку
равен уровню инфляции плюс 1-2%.
[11] Или числами ГКО, если
вспомнить 17 августа 1998 г.
[12] См. в конце книги
информацию о фирме СофтЛайн, официальном представителе MathSoft на
российском рынке (рекламная пауза!).
