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
|
// Copyright INRIA
//
// Van der Pol ODE
//
mode(-1)
xbasc()
titlepage(["Van der Pol:";...
"dy1/dt= y2";..
"dy2/dt=5*(1-y1^2)*y2-y1"]);
halt(); xbasc()
mode(1)
oldln=lines();
lines(0)
deff("yprim=f(t,y)",..
["yprim1=y(2)";..
"yprim2=mu*(1-y(1)^2)*y(2)-y(1)";..
"yprim=[yprim1;yprim2]";])
mu=5;
xmin=-3; xmax=3; ymin=-8; ymax=8;
xr=xmin:0.5:xmax; yr=ymin:1:ymax;
xbasc()
xset("pixmap",1)
fchamp(f,1,xr,yr)
xset("font",2,12)
xtitle([""],["y1"],["y2"])
xset("wshow")
xset("alufunction",6)
xset("thickness",2)
t0=0; dt=0.05; tmax=15;
t=t0:dt:tmax;
oldx0=10*xmax; oldy0=10*ymax;
dx=0.1; dy=0.1;
rtol=0.0001; atol=rtol;
while (%t)
[b,x0,y0]=xclick();
if b==2 then break end;
if b==0 & x0>=xmin & x0<=xmax & y0>=ymin & y0<=ymax then
sol=ode([x0;y0],t0,t,rtol,atol,f);
plot2d(sol(1,:),sol(2,:),5,"000")
xset("wshow")
rep=[x0,y0,1];
while rep(3)<>0 then
rep=xgetmouse(0);
x0=rep(1); y0=rep(2);
if (xmin<x0 & x0<xmax & ymin<y0 & y0<ymax) &..
(abs(x0-oldx0)>=dx | abs(y0-oldy0)>=dy) then
plot2d(sol(1,:),sol(2,:),5,"000")
sol=ode([x0;y0],t0,t,rtol,atol,f);
plot2d(sol(1,:),sol(2,:),5,"000")
xset("wshow")
oldx0=x0; oldy0=y0;
end
end
xset("alufunction",1)
plot2d(sol(1,:),sol(2,:),5,"000")
xset("wshow")
xset("alufunction",6)
end
end
xset("pixmap",0)
xset("default")
lines(oldln(1))
|
