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

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

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

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. [info]stepanishchev пишет:

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

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

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

  2. [info]soonts пишет:

    Уравнения безусловно решаем:
    http://soonts.freehostia.com/Runge-Kutta.html
    Надеюсь, бесплатный opera browser ты не отнесёшь к “коробочным продуктам”.