← Back to branch summary

~ubuntu-branches/ubuntu/karmic/scilab/karmic

« back to all changes in this revision

Viewing changes to examples/misc-examples/ode1.sce

  • 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:
examples/misc-examples/ode1.sce
 
1
//Copyright INRIA
 
2
Eps=1.e-2
 
3
//     Example of use of ode function:
 
4
//     System to solve:
 
5
//     dy1/dt = -0.04*y1 + 1.e4*y2*y3
 
6
//     dy2/dt = 0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
 
7
//     dy3/dt = 3.e7*y2**2
 
8
// on the interval from t = 0.0 to t = 4.e10, with initial conditions
 
9
// y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
 
10
//     Scilab definition of rhs function:
 
11
deff('[ydot]=f(t,y)',...
 
12
     'f1=-0.04*y(1)+1e4*y(2)*y(3),...
 
13
      f3=3e7*y(2)*y(2),...
 
14
      ydot=[f1;-f1-f3;f3]')
 
15
//     Definition of jacobian
 
16
deff('[jac]=j(t,y)',...
 
17
     'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
 
18
      jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
 
19
      jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);')
 
20
//    y0= Initial state  t1= times at which sol. is computed
 
21
y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
 
22
//    Reference solution is:
 
23
yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 
24
//
 
25
//  1. fortran program wfex and wjex dynamically called
 
26
files=G_make(["/tmp/wfex.o"],"wfex.dll");
 
27
iwfex=link(files,"wfex");
 
28
files=G_make(["/tmp/wjex.o"],"wjex.dll");
 
29
iwjex=link(files,"wjex");
 
30
//     jacobian is not given
 
31
y1=ode(y0,t0,t1,'wfex');
 
32
 
 
33
//  2. fortran called type given (stiff), no jacobian
 
34
y2=ode('stiff',y0,t0,t1,'wfex');
 
35
//  3. fortran called type (stiff) given
 
36
y3=ode('stiff',y0,t0,t1,'wfex','wjex');
 
37
//   hot restart
 
38
[z,w,iw]=ode('stiff',y0,0,0.4,'wfex','wjex');
 
39
z=ode('stiff',z,0.4,4,'wfex','wjex',w,iw);
 
40
if max(z-y3(:,2))  > Eps then pause,end
 
41
 
 
42
[y1,w,iw]=ode(y0,t0,t1(1),'wfex');
 
43
y2=ode(y0,t1(1),t1(2:nt),'wfex',w,iw);
 
44
if max([y1 y2]-yref)  > Eps then pause,end
 
45
 
 
46
[y1,w,iw]=ode(y0,t0,t1(1),'wfex','wjex');
 
47
y2=ode(y0,t1(1),t1(2:nt),'wfex','wjex',w,iw);
 
48
if max([y1 y2]-yref)  > Eps then pause,end
 
49
 
 
50
// Unlink wfex and wjex
 
51
 
 
52
ulink(iwfex);
 
53
ulink(iwjex);
 
54
//   Setting tolerances
 
55
atol=[0.001,0.0001,0.001];rtol=atol;
 
56
//   
 
57
//  4. type given , scilab lhs ,jacobian not passed
 
58
y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
 
59
if maxi(y4(:,1)-yref(:,1))  > Eps then pause,end
 
60
 
 
61
//  5. type non given,  scilab rhs and jacobian functions
 
62
y5=ode(y0,t0,t1,f,j);
 
63
if max(y5-yref)  > Eps then pause,end
 
64
//  6. type given (stiff),rhs and jacobian  by scilab
 
65
y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
 
66
if max(y6-yref)  > Eps then pause,end
 
67
 
 
68
//  7. matrix rhs, type given (adams),jacobian not passed
 
69
// 
 
70
a=rand(3,3);ea=expm(a);
 
71
deff('[ydot]=f(t,y)','ydot=a*y')
 
72
t1=1;y=ode('adams',eye(a),t0,t1,f);
 
73
if max(ea-y)  > Eps then pause,end
 
74
 
 
75
 
 
76
 
 
77