# Error drawing curve within slopefield

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;

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);
``````

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;

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:

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

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 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;
}
if(dx < 0) {
factor=(step+dx)/step;
g=g..{S3}cp0{S3};
break;
}
if(dy < 0) {
factor=(step+dy)/step;
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;

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