☰

坡度场内误差绘制曲线

tag-icon • tag-icon asymptote
坡度场内误差绘制曲线

请注意以下几行 Asymptote 代码中通过点 B 的绿色曲线和通过点 C 的红色曲线在 xmin 附近的图形中的错误。有没有办法纠正这个问题？

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

在此处输入图片描述

答案1

这是增量构建路径的一个已知问题（请guide参阅Asymptote 用户指南第 6.2 节)。我找不到忘记路径控制点的预定义方法，但你总是可以实现一个函数，该函数将从旧路径点构建新路径。convpath()以下代码中的函数就是这样做的：

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

结果是：

在此处输入图片描述

答案2

只是给出一个补充和不同的解决方案。我认为这个问题与方法和数值问题有关。

该例程curve使用 Runge-Kutta 四阶法来构造曲线和控制点（或多或少是每个点的导数）。由于是y=-2两个不同区域之间的边界，如果步长不够小，则S3计算的（或 RK4 方法中的 k4）位于减小的区域，而计算的点位于增大的区域，因此会出现这种奇怪的行为。您还可以观察到步长较小的明显差异。

我curve使用可选参数修改了例程stepfraction，请参阅mycurve下面代码中的例程。

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

结果

在此处输入图片描述

相关内容