← Back to branch summary

~ubuntu-branches/ubuntu/hoary/scilab/hoary

« back to all changes in this revision

Viewing changes to demos/ode/ode2.dem

  • Committer: Bazaar Package Importer
  • Author(s): Torsten Werner
  • Date: 2002-03-21 16:57:43 UTC
  • Revision ID: james.westby@ubuntu.com-20020321165743-e9mv12c1tb1plztg
Tags: upstream-2.6
Import upstream version 2.6

Show diffs side-by-side

added added

removed removed

Lines of Context:
demos/ode/ode2.dem
 
1
// Copyright INRIA
 
2
//chemical process (stiff)
 
3
mode(-1)
 
4
xbasc();
 
5
titlepage(['Integration of ODE:';...
 
6
        'dy1/dt=-0.04*y1 + 1d4*y2*y3';...
 
7
        'dy3/dt= 3d7*y2*y2';...
 
8
        'dy2/dt= -dy1/dt - dy3/dt';...
 
9
        'finding points such that';...
 
10
        'y1=1.e-4 or y3=1.e-2'])
 
11
mode(1)
 
12
// Equations definition
 
13
deff('[yd]=chem(t,y)',[
 
14
    'yd(1)=-0.04*y(1) + 1d4*y(2)*y(3);';
 
15
    'yd(3)= 3d7*y(2)*y(2);';
 
16
    'yd(2)= -yd(1) - yd(3);'])
 
17
 
 
18
// Integration
 
19
t=[1.d-5:0.02:.4 0.41:.1:4 40 400 4000 40000 4d5 4d6 4d7 4d8 4d9 4d10];
 
20
rtol=1.d-4;atol=[1.d-6;1.d-10;1.d-6];
 
21
y=ode([1;0;0],0,t,rtol,atol,chem);
 
22
// Visualisation
 
23
halt();xbasc();
 
24
rect=[1.d-5,-0.1,1d11,1.1];
 
25
plot2d1("oln",t',(diag([1 10000 1])*y)',(1:3),"111",' y1@10000 y2@y3',rect)
 
26
halt();
 
27
// Add surface condition
 
28
nt=prod(size(t));
 
29
deff('[y]=surf(t,x)','y=[x(1)-1.e-4;x(3)-1.e-2]')
 
30
 
 
31
// First root
 
32
[y,rd,w,iw]=ode('root',[1;0;0],0,t,rtol,atol,chem,2,surf);rd
 
33
 
 
34
while rd<>[] then 
 
35
  [nw,ny]=size(y);
 
36
  k=find(rd(1)>t(1:nt-1)&rd(1)<t(2:nt));
 
37
  // Visualisation       
 
38
  write(%io(2),[rd(1);y(:,ny)]','(''t='',e10.3,'' y='',3(e10.3,'',''))')
 
39
  plot2d1("oln",rd(1)',(diag([1 10000 1])*y(:,ny))',[-3,-3,-3],"000");
 
40
  // Next root
 
41
  [y,rd,w,iw]=ode('root',[1;0;0],rd(1),t(k+1:nt),rtol,atol,chem,2,surf,w,iw);
 
42
end
 
43