Предлагаю сравнить несколько языков программирования (в широком смысле) по критерию их удобства для разнообразных расчетов. Даю вводную: из коробки извлекается компьютер (или вообще какое-нибудь вычислительное устройство, желательно с экраном), и необходимо в кратчайший срок произвести на нем сравнительно несложные вычисления.
Чтобы не извращаться, возьмем возникающее из физики уравнение Ван-дер-Поля. Речь идет о «модификации» уравнения осциллятора
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):
Предлагаю всем желающим попробовать построить такой же график любыми доступными с компьютером «из коробки» средствами. К таковым, например, можно причислить Excel — офисный пакет довольно часто является предустановленным, Matlab же к «коробочным» не относится. Разнообразные варианты Basic — от Sinclair Basic до VBScript, Javascript, различные варианты интерпретатора командной строки — только приветствуются. Интересно было бы взглянуть на программу для какого-нибудь программируемого калькулятора (при отсутствии графического экрана достаточно просто «распечатки» значений). Сомневаюсь, что «старинные» МК-52 или МК-61 «потянут» это, памяти у них маловато (хотя метод Ньютона для них реализовывали), а вот для МК-152 такая задача — в самый раз.
PS Я свой график строил на QBasic из комплекта MS-DOS 6.22
![[info]](http://stat.livejournal.com/img/userinfo.gif){kind=small}
stepanishchev, естественно, с использованием МК-152. Время составления программы — порядка получаса, время счета — 2 минуты 15 секунд. График бодро выводится по мере расчетов :)
