Фрактал «магнитный маятник» - Основные уравнения

Основные уравнения

Движение маятника вычисляется путем двойного интегрирования по ускорениям, действующим на маятник. Обычно бы речь шла не об ускорениях, а о силах. Согласно первому закону Ньютона, сила, необходимая для перемещения тела, равна произведению массы на ускорение. Это уравнение было решено для ускорения.

 

Так как начальные условия дают начальное положение и начальную скорость (считающуюся нулевой), надо лишь вычислить ускорения. Для простоты масса считается равной одной единице массы. Как правило, моделирование особо не заботится о физических единицах. Это не проблема, поскольку использование реальных единиц лишь наложило бы коэффициенты пересчета на параметры. Следующие уравнения перечисляют все ускорения, важные для моделирования:

Код для прослеживания траектории маятника

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

  пока прослеживается траектория маятника
    положение += скорость
    ускорение = 0

    для всех источников силы
    ускорение += ускорение, вызванное источником

    если (маятник рядом с источником, и скорость маленькая) то
      магнит останова = индекс источника
      прервать
    конец ветвления if
    конец цикла for

    ускорение -= ускорение, вызванное трением
    скорость     += ускорение
    длина траектории += длина(скорость)

    сохранить магнит останова
    сохранить длину траектории
  конец цикла while

Реализация кода для прослеживания траектории, выбранной маятником, на C++ выглядит так:

for (int ct=0; ct<m_nMaxSteps && bRunning; ++ct)
{
  // вычислить новое положение
  pos[0] += vel[0]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[0] - 1.0 / 6.0 * (*acc_p)[0]);
  pos[1] += vel[1]*dt + sqr(dt) * (2.0/3.0 * 
           (*acc)[1] - 1.0 / 6.0 * (*acc_p)[1]);

  (*acc_n) = 0.0;  // сбросить ускорение

  // Вычислить силу, речь идет о силах, пропорциональных
  // расстоянию или обратному квадрату расстояния
  for (std::size_t i=0; i<src_num; ++i)
  {
    const Source &src( m_vSources[i] );
    r = pos - src.pos;
    if (src.type==Source::EType::tpLIN)
    {
      //---------------------------------------
      // Закон Гука:          _
      //         _             r         _   
      //  m * a = - k * |r| * --- = -k * r
      //                      |r|
      //
      (*acc_n)[0] -= src.mult * r[0];
      (*acc_n)[1] -= src.mult * r[1];
    } 
    else 
    {     
      //---------------------------------------
      // Магнитные силы: _
      //      _         r
      //  m * a = k * -----
      //               |r?|
      //
      double dist( sqrt( sqr(src.pos[0] - pos[0]) +  
                         sqr(src.pos[1] - pos[1]) + 
                         sqr(m_fHeight) ) );
      (*acc_n)[0] -= (src.mult / (dist*dist*dist)) * r[0]; 
      (*acc_n)[1] -= (src.mult / (dist*dist*dist)) * r[1]; 
    }

    // Проверить условие прерывания
    if (ct>m_nMinSteps && abs(r)<src.size && abs(vel)<m_fAbortVel)
    {
      bRunning = false;
      stop_mag = (int)i;
      break;
    }
  } // для всех источников

  //--------------------------------------------------------------
  // 3.) Трение пропорционально скорости 
  (*acc_n)[0] -= vel[0] * m_fFriction;
  (*acc_n)[1] -= vel[1] * m_fFriction;    

  //--------------------------------------------------------------
  // 4.) Интегрирование Бимена для скоростей
  vel[0] += dt * ( 1.0/3.0 * (*acc_n)[0] + 5.0/6.0 * 
                 (*acc)[0] - 1.0/6.0 * (*acc_p)[0] );
  vel[1] += dt * ( 1.0/3.0 * (*acc_n)[1] + 5.0/6.0 * 
                 (*acc)[1] - 1.0/6.0 * (*acc_p)[1] );

  //--------------------------------------------------------------
  // 5.) перевернуть буфер ускорений
  tmp = acc_p;
  acc_p = acc;
  acc   = acc_n;
  acc_n = tmp;
};