Error drawing curve within slopefield

It's a known issue with paths which are built incrementally (see the guide item of section 6.2 of the Asymptote user guide). I can't find a predefined way of forgetting the path control points, but you can always implement a function which would construct a new path from the old path points. The convpath() function in the following code does just that:

import graph;
import slopefield;
import fontsize;

defaultpen(fontsize(9pt));
size(300);
real dy(real x,real y) {return (y-1)^2*(y+2);}
real xmin=-2, xmax=5;
real ymin=-4, ymax=3;

add(slopefield(dy,(xmin,ymin),(xmax,ymax),21,black+0.5bp,Arrow));
draw((-2,1)--(5,1), blue+1bp);
draw((-2,-2)--(5,-2), blue+1bp);
pair B=(-1.0,-2.5);
pair C=(0.0,-1.5);
pair D=(0.0,2.5);

path convpath(path p) {
    guide res;
    for(int t=0; t<=length(p); ++t)
        res = res .. point(p, t);
    return res;
}

draw(convpath(curve(B,dy,(xmin,ymin),(xmax,ymax))),green+1bp);
draw(convpath(curve(C,dy,(xmin,ymin),(xmax,ymax))),red+1bp);
draw(curve(D,dy,(xmin,ymin),(xmax,ymax)),blue+1bp);

xaxis(YEquals(ymin),xmin,xmax,LeftTicks());
xaxis(YEquals(ymax),xmin,xmax);
yaxis(XEquals(xmin),ymin,ymax,RightTicks());
yaxis(XEquals(xmax),ymin,ymax);

The result is:

enter image description here


Just to give a complement and a different solution. I think that the problem is related to the method and a numerical issue.

The routine curve uses Runge-Kutta 4th order method to construct the curve and also the control point (more or less it is the derivative in each point). Since y=-2 is the frontier between two different zones, if the stepsize is not sufficiently small, the S3 (or k4 in the RK4 method) computed is in the decreasing area while the computed point is in the increasing area and as a consequence you have this strange behavior. You can also observe a visible difference with a smaller stepsize.

I modify the curve routine with an optionnal stepfraction argument, see mycurve routine in the following code.

    import graph;
    import slopefield;
    import fontsize;


    path mycurve(pair c, real f(real,real), pair a, pair b, real stepfr=0.05) 
    {
      real step=stepfr*(b.x-a.x);     
      real halfstep=0.5*step;
      real sixthstep=step/6;

      path follow(real sign) {
        pair cp=c;
        guide g=cp;
        real dx,dy;
        real factor=1;
        do {
          real slope;
          pair S(pair z) {
            slope=f(z.x,z.y);
            return factor*sign/sqrt(1+slope^2)*(1,slope);
          }
          pair S3;
          pair advance() {
            pair S0=S(cp);
            pair S1=S(cp+halfstep*S0);
            pair S2=S(cp+halfstep*S1);
            S3=S(cp+step*S2);
            pair cp0=cp+sixthstep*(S0+2S1+2S2+S3);
            dx=min(cp0.x-a.x,b.x-cp0.x);
            dy=min(cp0.y-a.y,b.y-cp0.y);
            return cp0;
          }
          pair cp0=advance();
          if(dx < 0) {
            factor=(step+dx)/step;
            cp0=advance();
            g=g..{S3}cp0{S3};
            break;
          }
          if(dy < 0) {
            factor=(step+dy)/step;
            cp0=advance();
            g=g..{S3}cp0{S3};
            break;
          }
          cp=cp0;
          g=g..{S3}cp{S3};
        } while (dx > 0 && dy > 0);
        return g;
      }

      return reverse(follow(-1))& follow(1);
    }


    path mycurve(pair c, real f(real), pair a, pair b, real stepfr=0.05)
    {
      return mycurve(c,new real(real x, real y){return f(x);},a,b,stepfr);
    }


    path convpath(path p) {
        guide res;
        for(int t=0; t<=length(p); ++t)
            res = res .. point(p, t);
        return res;
    }

    defaultpen(fontsize(9pt));
    size(300);
    real dy(real x,real y) {return (y-1)^2*(y+2);}
    real xmin=-2, xmax=5;
    real ymin=-4, ymax=3;

    add(slopefield(dy,(xmin,ymin),(xmax,ymax),21,black+0.5bp,Arrow));
    draw((-2,1)--(5,1), blue+1bp);
    draw((-2,-2)--(5,-2), blue+1bp);
    pair B=(-1.0,-2.5);
    pair C=(0.0,-1.5);
    pair D=(0.0,2.5);

    draw(convpath(curve(B,dy,(xmin,ymin),(xmax,ymax))),green);
    draw(mycurve(B,dy,(xmin,ymin),(xmax,ymax),0.02),lightgreen);
    draw(mycurve(C,dy,(xmin,ymin),(xmax,ymax),0.02),lightred);
    draw(convpath(curve(C,dy,(xmin,ymin),(xmax,ymax))),red);
    draw(mycurve(D,dy,(xmin,ymin),(xmax,ymax)),blue+1bp);


    xaxis(YEquals(ymin),xmin,xmax,LeftTicks());
    xaxis(YEquals(ymax),xmin,xmax);
    yaxis(XEquals(xmin),ymin,ymax,RightTicks());
    yaxis(XEquals(xmax),ymin,ymax);  

And the result

enter image description here

Tags:

Asymptote