Re: Isoclines

["Andrew D. Hwang" <ahwang@mathcs.holycross.edu>]

On Mon, 6 Jun 2005, Gunnar wrote:

> Is it possible to draw isoclines with ePiX in some easy way?
>
> Just to be clear, with isoclines to an ODE y'=f(x,y) I mean the curves on
> which y' has a constant value C, and the isoclines are drawn for different
> values of C.
>
If I understand correctly, you want to draw level curves of f(x,y).  One
trick that's worked for me is to compute J(grad f), the gradient of f
rotated a quarter turn, then to draw solution curves of the resulting
planar ODE. It's not elegant; you'll have to choose the initial points by
hand, and may have to adjust the plot intervals manually to avoid
"wrapping", but depending how complicated f is, the results should look
reasonably good.

A short sample is attached.

Regards,
Andy

Andrew D. Hwang			ahwang@mathcs.holycross.edu
Department of Math and CS	http://math.holycross.edu/~ahwang
College of the Holy Cross	(508) 793-2458 (Office: 320 Swords)
Worcester, MA, 01610-2395	(508) 793-3530 (fax)

/* -*-ePiX-*- */
#include "epix.h"
using namespace ePiX;

int arms=6;

P f(double x, double y)
{
  return P(1, x*(x*x - 3*y*y));
}

// J(grad f)
P F(double x, double y)
{
  return P(6*x*y, 3*(x*x-y*y));
}

void isocline(P arg)
{
  // kludge to plot forward and backward; ode_plot seems to be buggy
  ode_plot(F, arg, 0,  1.5, 120);
  ode_plot(F, arg, 0, -1.5, 120);
}

int main() 
{
  bounding_box(P(-2,-2),P(2,2));
  unitlength("1cm");
  picture(10,10);

  begin();

  crop();
  slope_field(f, P(x_min, y_min), P(x_max, y_max), 16, 16);

  degrees();
  bold();
  red();

  for (int i=1; i < 4; ++i)
    for (int j=0; j< arms; ++j)
      isocline(polar(0.5*i, j*360.0/arms));

  end();
}