1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
import graph_settings;
real stepfraction=0.05;
picture slopefield(real f(real,real), pair a, pair b,
int nx=nmesh, int ny=nx,
real tickfactor=0.5, pen p=currentpen)
{
picture pic;
real dx=(b.x-a.x)/nx;
real dy=(b.y-a.y)/ny;
real step=0.5*tickfactor*min(dx,dy);
for(int i=0; i <= nx; ++i) {
real x=a.x+i*dx;
for(int j=0; j <= ny; ++j) {
pair cp=(x,a.y+j*dy);
real slope=f(cp.x,cp.y);
real mp=step/sqrt(1+slope^2);
draw(pic,(cp.x-mp,cp.y-mp*slope)--(cp.x+mp,cp.y+mp*slope),p);
}
}
return pic;
}
picture slopefield(real f(real), pair a, pair b,
int nx=nmesh, int ny=nx, pen p=currentpen)
{
return slopefield(new real(real x, real y) {return f(x);},a,b,nx,ny,p);
}
path curve(pair c, real f(real,real), pair a, pair b)
{
real step=stepfraction*(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 curve(pair c, real f(real), pair a, pair b)
{
return curve(c,new real(real x, real y){return f(x);},a,b);
}
|
