7

Note the error in the graphs of the green curve passing through point B and the red curve passing through point C near xmin for the following lines of Asymptote code. Is there a way to correct for this?

    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);
    draw(curve(B,dy,(xmin,ymin),(xmax,ymax)),green+1bp);
    draw(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);  

enter image description here

7

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

  • Thank you! New to using Asymptote, I suspected that this was a known issue, but didn't know where to look. Very good. – Stephen Clark Dec 25 '17 at 22:33
3

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.