Skip to content

Single Systems

Stephen Berry edited this page Jun 21, 2018 · 5 revisions

There are a number of ways to set up a system of ordinary differential equations in Ascent.

Here are examples of how to set up and solve a Lorenz attractor via a functor, function, and lambda function.

Functor

struct Lorenz
{
   void operator()(const state_t& x, state_t& xd, const double)
   {
      static constexpr double sigma = 10.0;
      static constexpr double R = 28.0;
      static constexpr double b = 8.0 / 3.0;

      xd[0] = sigma * (x[1] - x[0]);
      xd[1] = R * x[0] - x[1] - x[0] * x[2];
      xd[2] = -b * x[2] + x[0] * x[1];
   }
};

int main()
{
   state_t x = { 10.0, 1.0, 1.0 };
   double t = 0.0;
   double dt = 0.01;
   double t_end = 10.0;

   RK4 integrator;
   Recorder recorder;

   while (t < t_end)
   {
      recorder({ t, x[0], x[1], x[2] });
      integrator(Lorenz(), x, t, dt);
   }

   recorder.csv("lorenz", { "t", "x0", "x1", "x2" });

   return 0;
}

Function

void lorenz(const state_t& x, state_t& xd, const double)
{
   static constexpr double sigma = 10.0;
   static constexpr double R = 28.0;
   static constexpr double b = 8.0 / 3.0;

   xd[0] = sigma * (x[1] - x[0]);
   xd[1] = R * x[0] - x[1] - x[0] * x[2];
   xd[2] = -b * x[2] + x[0] * x[1];
}

int main()
{
   state_t x = { 10.0, 1.0, 1.0 };
   double t = 0.0;
   double dt = 0.01;
   double t_end = 10.0;

   RK4 integrator;
   Recorder recorder;

   while (t < t_end)
   {
      recorder({ t, x[0], x[1], x[2] });
      integrator(lorenz, x, t, dt);
   }

   recorder.csv("lorenz", { "t", "x0", "x1", "x2" });

   return 0;
}

Lambda Function

int main()
{
   auto lorenz = [](const state_t& x, state_t& xd, const double)
   {
      static constexpr double sigma = 10.0;
      static constexpr double R = 28.0;
      static constexpr double b = 8.0 / 3.0;

      xd[0] = sigma * (x[1] - x[0]);
      xd[1] = R * x[0] - x[1] - x[0] * x[2];
      xd[2] = -b * x[2] + x[0] * x[1];
   };

   state_t x = { 10.0, 1.0, 1.0 };
   double t = 0.0;
   double dt = 0.01;
   double t_end = 10.0;

   RK4 integrator;
   Recorder recorder;

   while (t < t_end)
   {
      recorder({ t, x[0], x[1], x[2] });
      integrator(lorenz, x, t, dt);
   }

   recorder.csv("lorenz", { "t", "x0", "x1", "x2" });

   return 0;
}

3D Plot of Results:

Clone this wiki locally