Решаем уравнение бурбулятора

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

Чтобы не извращаться, возьмем возникающее из физики уравнение Ван-дер-Поля. Речь идет о «модификации» уравнения осциллятора

y» + ay’ + y = 0.

Как легко показать (ну-ка, кто умеет решать дифференциальные уравнения?) при a>0 колебания будут затухающими, а при a<0 - неустойчивыми, то есть система пойдет вразнос. При a=0 получится малоинтересное уравнение идеального колебательного контура, которого в природе не бывает. Параметр a играет здесь роль сопротивления. Теперь подставим вместо константы a изменяющуюся величину (например, включив в цепь триод). В качестве простого "идеализированного" триода выберем a = y2 — 1. Тогда уравнение превратится в

y» + (y2-1)y’ + y = 0,

или, после стандартного понижения порядка — в систему двух уравнений:

y’1=y2
y’2=(1 — y12)y2 — y1

Как показал в 1920-1926 году Ван-дер-Поль, у такой системы сущеcтвует устойчивое периодическое решение, к которому сходятся все остальные решения. Это не может не радовать, так как система возникла из математического описания генератора незатухающих колебаний на триоде (кстати, для 20-х — вполне себе чудо техники).

Необходимые же вычисления заключаются в том, чтобы численно промоделировать это уравнение. Предлагаю воспользоваться методом Рунге-Кутта в его простейшей форме. Пусть у нас есть система дифференциальных уравнений

y’1=f1(x, y1, y2)
y’2=f2(x, y1, y2),

мы умеем вычислять функции f1 и f2 в любой точке, и нам известны значения y1, y2 в точке x0. Тогда мы можем приближенно вычислить их значения в точке x1 = x0 + h, где h достаточно мало, по следующим формулам (в «Кратком физико-математическом справочнике» Аленицына, Бутикова и Кондратьева 1990 года издания они приведены, как «формулы Рунге-Кутта»):

p11 = hf1(x0, y1, y2)
p12 = hf2(x0, y1, y2)

p21 = hf1(x0 + h/2, y1 + p11/2, y2 + p12/2)
p22 = hf2(x0 + h/2, y1 + p11/2, y2 + p12/2)

p31 = hf1(x0 + h/2, y1 + p21/2, y2 + p22/2)
p32 = hf2(x0 + h/2, y1 + p21/2, y2 + p22/2)

p41 = hf1(x0 + h, y1 + p31, y2 + p32)
p42 = hf2(x0 + h, y1 + p31, y2 + p32)

x1 = x0 + h
y1(x1) ≈ y01 + (p11 + 2 p21 + 2 p31 + p41) / 6
y2(x1) ≈ y02 + (p12 + 2 p22 + 2 p32 + p42) / 6

Конечно, в современной вычислительной практике применяются гораздо более сложные варианты метода Рунге-Кутта, например, формулы Дормана-Принса, но для грубых расчетов хватит и этого. Не будем усложнять задачу такими «бонусами», как автоматический выбор длины шага, а ограничимся этими формулами. Для упрощения написания программы я «раскрыл» вектор (y1, y2). Для вычисления y1(x1) и y2(x1) достаточно лишь применить эти формулы в том порядке, в котором они выписаны. Вычисления надо повторять, пока не исчерпается отрезок, на котором нас интересуют значения y(x).

Теперь для расчетов требуются лишь начальные условия — известные в нулевой момент времени y1 и y2, а также величина шага и длина отрезка интегрирования. Предлагаю выбрать такие начальные условия:

x0 = 0,
y1(0) = 0,
y2(0) = 0,0001

Очень небольшое значение y2 — необходимый для запуска генератора «толчок». Я утверждаю, что при таких условиях y1(x) довольно быстро «стабилизируется» и колебания станут незатухающими, с постоянной частотой. Для того, чтобы убедиться в этом, предлагаю проделать вычисления с шагом 0,1 вплоть до x = 64. Получим довольно симпатичный график (масштаб по y — в десять раз меньше, чем по x):

graph

Предлагаю всем желающим попробовать построить такой же график любыми доступными с компьютером «из коробки» средствами. К таковым, например, можно причислить Excel — офисный пакет довольно часто является предустановленным, Matlab же к «коробочным» не относится. Разнообразные варианты Basic — от Sinclair Basic до VBScript, Javascript, различные варианты интерпретатора командной строки — только приветствуются. Интересно было бы взглянуть на программу для какого-нибудь программируемого калькулятора (при отсутствии графического экрана достаточно просто «распечатки» значений). Сомневаюсь, что «старинные» МК-52 или МК-61 «потянут» это, памяти у них маловато (хотя метод Ньютона для них реализовывали), а вот для МК-152 такая задача — в самый раз.

PS Я свой график строил на QBasic из комплекта MS-DOS 6.22

Решаем уравнение бурбулятора: 2 комментария

  1. Решение на ЭКВМ Электроника МК-152/МК-161 по ссылке:
    http://mk.semico.ru/dr_info25.htm

    Длина программа 248 байт вместе с выводом графики и записью всех точек в электронный блокнот. Время составления программы — порядка получаса. Время работы около минуты — с отрисовкой графики по ходу работы. Если выводить график один раз, то будет значительно быстрее.

    Программа не оптимизирована. Содержит несколько лишних операторов и блок инициализации начальных значений.

Добавить комментарий для ext:stepanishchev@lj Отменить ответ

Ваш адрес email не будет опубликован. Обязательные поля помечены *