За последние 40 лет применение компьютеров возродило старую и тяжелую задачу классической механики. Все началось с системы трех связанных ДУ, придуманных Эдом Лоренцом (Ed Lorentz) при попытках моделирования погоды.

• restart;with(plots):
• Eq1:=diff(x(t),t) = sigma*(y(t)-x(t));
• Eq2:=diff(y(t),t) = r*x(t)-y(t)-x(t)*z(t);
• Eq3:=diff(z(t),t) = x(t)*y(t)-b*z(t);

где константы σ, r, b можно выбирать как хочется. В большинстве случаев выбор параметров не давал ничего интересного, но Лоренц открыл, что есть диапазон параметров, в котором решения ведут себя очень странно. Если выберем

• sigma:=10;b:=8/3;

и будем менять r от 0 до ~30, происходит переход от скучного ординарного поведения к странному, которое революционизировало много областей науки, в том числе физику.

Для численного решения этого набора ДУ с особым выбором r применим:

• r:=10;s:=dsolve({Eq1,Eq2,Eq3,x(0)=1,y(0)=1,z(0)=1},
  {x(t),y(t),z(t)},type=numeric);

и с помощью odeplot посмотрим, как выглядит решение:

• odeplot(s,[[t,x(t)],[t,y(t)],[t,z(t)]],0..30,numpoints=600);

Запустим группу для значений r: 10, 20, 28, 35, 45. При r выше ≈ 25 колебания создают рисунок с не повторяющейся структурой. (Запустите odeplot на более долгое время, чтобы отсортировать колебания.) А вот что совсем плохо: если взять чуть другие начальные условия, получается совсем иная картинка (это знаменитый эффект бабочки. Чтобы его увидеть, возьмите r = 50 и сравните результаты при начальных условиях x(0) = 1.0 и x(0) = 1.0000001, обращая особое внимание на то, как выглядят зеленый и красный следы вблизи t = 30). Это совсем неожиданно. Управляющие этим поведением уравнения полностью определены и достаточно просты. Как может из простой вещи появиться случайный выброс? Это явно удивительно, и такое поведение называется динамическим хаосом.

Другой способ. Вместо графика x(t) рисуем с помощью команды графики spacecurve трехмерную траекторию кривой решения [x(t), y(t), z(t)]. Для этого делается цикл повторяющихся вызовов dsolve, создающих последовательность трехмерных графиков вдоль траектории системы. В spacecurve передается последовательность точек, как в примере ниже.

Зададим для r0 начальные значения (x, y, z), начальное время t1, конечное время t2 и число точек вдоль траектории numpoints, затем вычислим шаг по времени между двумя последовательными точками:

• r0:=[1,1,1];t1:=0;t2:=10;numpoints:=400;
  dt:=(t2-t1)/numpoints;

Зададим параметр системы r:

• r:=50;

Шаг системы включает повторяющиеся вызовы dsolve, каждый из которых своим начальным условием имеет конечное условие предыдущего шага. Это сделано циклом for-do-enddo (см. главу Программирование). Учтите, что вычисления длительные и придется ждать.

• for i from 1 to numpoints do
• s:=dsolve({Eq1,Eq2,Eq3,x(0)=r0[1],y(0)=r0[2],z(0)=r0[3]},
  {x(t),y(t),z(t)},type=numeric):
• t0:=t1+(i-1)*dt:
• tf:=t0+dt:
• sol:=s(tf-t0):

После каждого шага сохраним (x, y, z) в rf[i]:

• rf[i]:=[rhs(op(2,sol)),rhs(op(3,sol)),rhs(op(4,sol))]:

Получим следующее начальное условие:

• r0:=rf[i]:

Окончание цикла, возврат назад и повторное исполнение:

• enddo:

Поставим rf в последовательность, чтобы использовать со spacecurve:

• L:=[seq(rf[i],i=1..numpoints)]:

Построим траекторию:

• spacecurve(L,orientation=[45,60]);

Можно поставить курсор и покрутить график или анимировать вращение в цикле.

Сделаем 100 пространственных кривых, каждую с немного другим углом обзора. Затем анимируем эту последовательность так, что траектория поворачивается на экране, облегчая ее рассмотрение в 3D-:

• for i from 1 to 100 do
• theta :=45:phi:=i*3.6:
• p[i]:=spacecurve(L,orientation=[theta,phi]):
• od:

Нарисуем анимацию.

Когда появится картинка, щелкните по ней, чтобы возникли элементы управления анимацией на панели управления:

• plots[display]([seq(p[i],i=1..100)],insequence=true);

Конечно, можно щелкнуть на одной кривой и использовать курсор для ее вращения, но лучше посмотреть анимацию пространственной фигуры.