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
|
// Copyright INRIA
//chemical process (stiff)
mode(-1)
xbasc();
titlepage(['Integration of ODE:';...
'dy1/dt=-0.04*y1 + 1d4*y2*y3';...
'dy3/dt= 3d7*y2*y2';...
'dy2/dt= -dy1/dt - dy3/dt';...
'finding points such that';...
'y1=1.e-4 or y3=1.e-2'])
mode(1)
// Equations definition
deff('[yd]=chem(t,y)',[
'yd(1)=-0.04*y(1) + 1d4*y(2)*y(3);';
'yd(3)= 3d7*y(2)*y(2);';
'yd(2)= -yd(1) - yd(3);'])
// Integration
t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10];
rtol=1.d-4;atol=[1.d-6;1.d-10;1.d-6];
y=ode([1;0;0],0,t,rtol,atol,chem);
// Visualisation
halt();xbasc();
rect=[1.d-5,-0.1,1d11,1.1];
plot2d1("oln",t',(diag([1 10000 1])*y)',(1:3),"111",' y1@10000 y2@y3',rect)
halt();
// Add surface condition
nt=prod(size(t));
deff('[y]=surf(t,x)','y=[x(1)-1.e-4;x(3)-1.e-2]')
// First root
[y,rd,w,iw]=ode('root',[1;0;0],0,t,rtol,atol,chem,2,surf);rd
while rd<>[] then
[nw,ny]=size(y);
k=find(rd(1)>t(1:nt-1)&rd(1)<t(2:nt));
// Visualisation
write(%io(2),[rd(1);y(:,ny)]','(''t='',e10.3,'' y='',3(e10.3,'',''))')
plot2d1("oln",rd(1)',(diag([1 10000 1])*y(:,ny))',[-3,-3,-3],"000");
// Next root
[y,rd,w,iw]=ode('root',[1;0;0],rd(1),t(k+1:nt),rtol,atol,chem,2,surf,w,iw);
end
|
