← Back to branch summary

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

« back to all changes in this revision

Viewing changes to tests/matode.dia.ref

  • 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:
tests/matode.dia.ref
 
1
 
 
2
// Copyright INRIA
 
3
 
 
4
Leps=1.e-6;
 
5
 
 
6
//     dy1/dt = -0.04*y1 + 1.e4*y2*y3
 
7
 
 
8
//     dy2/dt =0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
 
9
 
 
10
//     dy3/dt = 3.e7*y2**2
 
11
 
 
12
// on the interval from t = 0.0 to t = 4.e10, with initial conditions
 
13
 
 
14
// y1 = 1.0, y2 = y3 = 0.  the problem is stiff.
 
15
 
 
16
//     definition of rhs
 
17
 
 
18
deff('[ydot]=f(t,y)',...
 
19
   'f1=-0.04*y(1)+1e4*y(2)*y(3),...
 
20
    f3=3e7*y(2)*y(2),...
 
21
    ydot=[f1;-f1-f3;f3]','n')
 
22
 
 
23
//     jacobian
 
24
 
 
25
deff('[jac]=j(t,y)',...
 
26
'jac(1,1)=-0.04;jac(1,2)=1.e4*y(3);jac(1,3)=1.e4*y(2),...
 
27
 jac(3,1)=0;jac(3,2)=6.e7*y(2);jac(3,3)=0;...
 
28
 jac(2,1)=-jac(1,1);jac(2,2)=-jac(1,2)-jac(3,2);jac(2,3)=-jac(1,3);','n')
 
29
 
 
30
//
 
31
 
 
32
y0=[1;0;0];t0=0;t1=[0.4,4];nt=size(t1,'*');
 
33
 
 
34
//    solution
 
35
 
 
36
yref=[0.9851721 0.9055180;0.0000339 0.0000224;0.0147940 0.0944596];
 
37
 
 
38
//
 
39
 
 
40
//  1. fortran called by fydot, without jacobian
 
41
 
 
42
y1=ode(y0,t0,t1,'fex');
 
43
 
 
44
if maxi(y1-yref) > Leps then bugmes();quit;end
 
45
 
 
46
//  2. fortran called by fydot, type given (stiff), no jacobian
 
47
 
 
48
y2=ode('stiff',y0,t0,t1,'fex');
 
49
 
 
50
if maxi(y2-yref) > Leps then bugmes();quit;end
 
51
 
 
52
//  3. fortran called by fydot , fjac, type given
 
53
 
 
54
y3=ode('stiff',y0,t0,t1,'fex','jex');
 
55
 
 
56
if maxi(y3-yref) > Leps then bugmes();quit;end
 
57
 
 
58
//   hot restart
 
59
 
 
60
[z,w,iw]=ode('stiff',y0,0,0.4,'fex','jex');
 
61
 
 
62
z=ode('stiff',z,0.4,4,'fex','jex',w,iw);
 
63
 
 
64
if maxi(z-y3(:,2)) > %eps then bugmes();quit;end
 
65
 
 
66
 
 
67
[y1,w,iw]=ode(y0,t0,t1(1),'fex');
 
68
 
 
69
y2=ode(y0,t1(1),t1(2:nt),'fex',w,iw);
 
70
 
 
71
if maxi([y1 y2]-yref) > Leps then bugmes();quit;end
 
72
 
 
73
 
 
74
[y1,w,iw]=ode(y0,t0,t1(1),'fex','jex');
 
75
 
 
76
y2=ode(y0,t1(1),t1(2:nt),'fex','jex',w,iw);
 
77
 
 
78
if maxi([y1 y2]-yref) > Leps then bugmes();quit;end
 
79
 
 
80
 
 
81
//   variation of tolerances
 
82
 
 
83
atol=[0.001,0.0001,0.001];rtol=atol;
 
84
 
 
85
//    externals
 
86
 
 
87
//  4. type given , scilab lhs ,jacobian not passed
 
88
 
 
89
y4=ode('stiff',y0,t0,t1(1),atol,rtol,f);
 
90
 
 
91
if maxi(y4(:,1)-yref(:,1)) > 0.01 then bugmes();quit;end
 
92
 
 
93
//  5. type non given, rhs and scilab jacobian
 
94
 
 
95
y5=ode(y0,t0,t1,f,j);
 
96
 
 
97
if maxi(y5-yref) > Leps then bugmes();quit;end
 
98
 
 
99
//  6. type given (stiff),rhs and jacobian  by scilab
 
100
 
 
101
y6=ode('stiff',y0,t0,t1,0.00001,0.00001,f,j);
 
102
 
 
103
if (y6-yref) > 2*0.00001 then bugmes();quit;end
 
104
 
 
105
//  7. matrix rhs, type given(adams),jacobian not passed
 
106
 
 
107
//
 
108
 
 
109
a=rand(3,3);ea=expm(a);
 
110
 
 
111
deff('[ydot]=f(t,y)','ydot=a*y')
 
112
Warning :redefining function: f                       
 
113
 
 
114
 
 
115
t1=1;y=ode('adams',eye(a),t0,t1,f);
 
116
 
 
117
if maxi(ea-y) > Leps then bugmes();quit;end
 
118
 
 
119
//
 
120
 
 
121
//   DAE's
 
122
 
 
123
//     dy1/dt = -0.04*y1 + 1.e4*y2*y3
 
124
 
 
125
//     dy2/dt =0.04*y1 - 1.e4*y2*y3 - 3.e7*y2**2
 
126
 
 
127
//       0.   = y1 + y2 + y3 - 1.
 
128
 
 
129
//  y1(0) = 1.0, y2(0) = y3(0) = 0.
 
130
 
 
131
// scilab macros
 
132
 
 
133
deff('[r]=resid(t,y,s)','...
 
134
r(1)=-0.04*y(1)+1.d4*y(2)*y(3)-s(1),...
 
135
r(2)=0.04*y(1)-1.d4*y(2)*y(3)-3.d7*y(2)*y(2)-s(2),...
 
136
r(3)=y(1)+y(2)+y(3)-1')
 
137
 
 
138
deff('[p]=aplusp(t,y,p)','...
 
139
p(1,1)=p(1,1)+1,...
 
140
p(2,2)=p(2,2)+1')
 
141
 
 
142
deff('[p]=dgbydy(t,y,s)','[-0.04,1.d4*y(3),1.d4*y(2);...
 
143
0.04,-1.d4*y(3)-6.d7*y(2),-1.d4*y(2);...
 
144
1,1,1]')
 
145
 
 
146
//        %ODEOPTIONS tests
 
147
 
 
148
//
 
149
 
 
150
rand('seed',0);rand('normal');
 
151
 
 
152
nx=20;A=rand(nx,nx);A=A-4.5*eye();
 
153
 
 
154
deff('y=f(t,x)','y=A*x')
 
155
Warning :redefining function: f                       
 
156
 
 
157
 
 
158
deff('J=j(t,x)','J=A')
 
159
Warning :redefining function: j                       
 
160
 
 
161
 
 
162
x0=ones(nx,1);t0=0;t=[1,2,3,4,5];
 
163
 
 
164
nt=size(t,'*');
 
165
 
 
166
eAt=expm(A*t(nt));
 
167
 
 
168
 
 
169
//        Test itask=%ODEOPTIONS(1)
 
170
 
 
171
 
 
172
//itask=1  --->  usual call (t=vector of instants, solution at all t)
 
173
 
 
174
//========================
 
175
 
 
176
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
177
 
 
178
jacflag=2;ml=-1;mu=-1;
 
179
 
 
180
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
181
 
 
182
xf=ode(x0,t0,t,f);   //lsoda
 
183
 
 
184
if norm(xf(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 
185
 
 
186
xfj=ode(x0,t0,t,f,j);   //lsoda with jacobian
 
187
 Warning: Jacobian external is given, but 
 
188
 not used!,  see %ODEOPTIONS(6)
 
189
 
 
190
 
 
191
if norm(xfj(:,nt)-eAt*x0)>Leps then bugmes();quit;end
 
192
 
 
193
 
 
194
 
 
195
//itask=2;   --->  solution at mesh points  ---> t=tmax
 
196
 
 
197
//========================
 
198
 
 
199
%ODEOPTIONS(1)=2;tmax=t(5);
 
200
 
 
201
xft=ode(x0,t0,tmax,f);
 
202
 
 
203
[p,q]=size(xft);
 
204
 
 
205
xlast=xft(2:nx+1,q);
 
206
 
 
207
if xft(1,q)<tmax then bugmes();quit;end
 
208
 
 
209
if norm(xlast-expm(A*xft(1,q))*x0)>Leps then bugmes();quit;end
 
210
 
 
211
 
 
212
//itask=3;   ---> solution at first mesh point beyond t=tmax
 
213
 
 
214
%ODEOPTIONS(1)=3;
 
215
 
 
216
x3=ode(x0,t0,tmax,f);
 
217
 
 
218
if norm(x3(2:nx+1)-xlast,1)>Leps then bugmes();quit;end
 
219
 
 
220
 
 
221
//itask=4; test with %ODEOPTIONS(2)=tcrit
 
222
 
 
223
%ODEOPTIONS(1)=4; //---> computation at all t and t>tcrit are not called
 
224
 
 
225
tcrit=2.5;%ODEOPTIONS(2)=tcrit;
 
226
 
 
227
chk=0;
 
228
 
 
229
deff('y=fcrit(t,x)',['if t<=tcrit then'
 
230
                      ' y=A*x;'
 
231
                      'else'
 
232
                      ' y=A*x;chk=resume(1);end'])
 
233
 
 
234
x42=ode(x0,t0,t,fcrit);
 
235
 Warning: integration up to tcrit
 
236
 
 
237
 
 
238
if chk==1 then bugmes();quit;end
 
239
 
 
240
[p,q]=size(x42);
 
241
 
 
242
if norm(x42(:,q)-ode(x0,t0,tcrit,f),1)>Leps then bugmes();quit;end
 
243
 
 
244
 
 
245
//itask=5; test with %ODEOPTIONS(2)=tcrit
 
246
 
 
247
%ODEOPTIONS(1)=5;  //---> computation at mesh points and t>tcrit are not called
 
248
 
 
249
%ODEOPTIONS(6)=2;  // Estimated jacobian
 
250
 
 
251
chk=0;
 
252
 
 
253
x52=ode(x0,t0,2.3,fcrit);
 
254
 
 
255
if chk==1 then bugmes();quit;end
 
256
 
 
257
[p,q]=size(x52);
 
258
 
 
259
if x52(1,q)>tcrit then bugmes();quit;end
 
260
 
 
261
 
 
262
//test of %ODEOPTIONS(3:5)=[h0,hmax,hmin]
 
263
 
 
264
%ODEOPTIONS(1)=1;
 
265
 
 
266
h0=0.0;hmax=0.1;hmin=0.0001;
 
267
 
 
268
%ODEOPTIONS(3:5)=[h0,hmax,hmin];
 
269
 
 
270
x35=ode(x0,t0,t,f);
 
271
 
 
272
if norm(x35-xf,1)>Leps then bugmes();quit;end
 
273
 
 
274
 
 
275
//test of %ODEOPTIONS(6)=jacflag
 
276
 
 
277
%ODEOPTIONS(6)=1;//Jacobian given
 
278
 
 
279
%ODEOPTIONS(3:5)=[0 0 0];
 
280
 
 
281
x61=ode('st',x0,t0,t,f,j);   //with Jacobian
 
282
 
 
283
if norm(x61-xf,1)>10*Leps then bugmes();quit;end
 
284
 
 
285
 
 
286
 
 
287
%ODEOPTIONS(6)=0; // jacobian nor called nor estimated
 
288
 
 
289
x60=ode('st',x0,t0,t,f,j);   //Jacobian not used (warning)
 
290
 Warning: Jacobian external is given, but 
 
291
 not used!,  see %ODEOPTIONS(6)
 
292
 
 
293
 
 
294
x60=ode('st',x0,t0,t,f);    //Jacobian not used
 
295
 
 
296
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
297
 
 
298
 
 
299
 
 
300
%ODEOPTIONS(6)=1;//Jacobian estimated
 
301
 
 
302
x60=ode('st',x0,t0,t,f)  ;
 
303
 Warning: No Jacobian external given but 
 
304
 one is required by %ODEOPTIONS(6) value!
 
305
 
 
306
 
 
307
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
308
 
 
309
 
 
310
//test of %ODEOPTIONS(6)=jacflag   (adams)
 
311
 
 
312
%ODEOPTIONS(6)=1;//with given Jacobian
 
313
 
 
314
x60=ode('ad',x0,t0,t,f,j) ;
 
315
 
 
316
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
317
 
 
318
 
 
319
 
 
320
%ODEOPTIONS(6)=0;// jacobian nor called nor estimated
 
321
 
 
322
x60=ode('ad',x0,t0,t,f,j);   //Jacobian not used (warning)
 
323
 Warning: Jacobian external is given, but 
 
324
 not used!,  see %ODEOPTIONS(6)
 
325
 
 
326
 
 
327
x60=ode('ad',x0,t0,t,f);    //Jacobian not used
 
328
 
 
329
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
330
 
 
331
 
 
332
// test lsoda
 
333
 
 
334
%ODEOPTIONS(6)=2;// jacobian  estimated
 
335
 
 
336
x60=ode(x0,t0,t,f);
 
337
 
 
338
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
339
 
 
340
 
 
341
%ODEOPTIONS(6)=1;//Jacobian estimated
 
342
 
 
343
x60=ode(x0,t0,t,f);
 
344
 Warning: No Jacobian external given but 
 
345
 one is required by %ODEOPTIONS(6) value!
 
346
 
 
347
 
 
348
if norm(x60-x61,1)>10*Leps then bugmes();quit;end
 
349
 
 
350
 
 
351
 
 
352
//   Banded Jacobian
 
353
 
 
354
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
355
 
 
356
//provisional values as default
 
357
 
 
358
jacflag=2;ml=-1;mu=-1;
 
359
 
 
360
 
 
361
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
362
 
 
363
jacflag=2;%ODEOPTIONS(6)=jacflag;   //Banded Jacobian, given
 
364
 
 
365
nx=20;A=diag(-[1:nx]);x0=ones(nx,1);
 
366
 
 
367
t0=0;t=[1,2,3,4,5];
 
368
 
 
369
for k=1:nx-1, A(k,k+1)=1;end
 
370
 
 
371
for k=1:nx-2, A(k,k+2)=-2;end
 
372
 
 
373
for k=1:nx-1, A(k+1,k)=-1;end
 
374
 
 
375
for k=1:nx-2, A(k+2,k)=2;end
 
376
 
 
377
for k=1:nx-3, A(k+3,k)=-3;end
 
378
 
 
379
deff('xd=f(t,x)','xd=A*x')
 
380
Warning :redefining function: f                       
 
381
 
 
382
 
 
383
ml=3;mu=2;
 
384
 
 
385
%ODEOPTIONS(11:12)=[ml,mu];
 
386
 
 
387
for i=1:nx;
 
388
    for j=1:nx;
 
389
if A(i,j)<>0 then J(i-j+mu+1,j)=A(i,j);end
 
390
end;end;
 
391
Warning :redefining function: j                       
 
392
 
 
393
 
 
394
// J is a ml+mu+1 x ny matrix.
 
395
 
 
396
// Column 1 of J is made of mu zeros followed by df1/dx1, df2/dx1, df3/dx1,
 
397
 
 
398
// i.e. 1 + ml possibly nonzeros entries.
 
399
 
 
400
// Column 2 of J is made of mu-1 zeros followed by df1/dx2, df2/dx2,0...
 
401
 
 
402
// etc...
 
403
 
 
404
%ODEOPTIONS(6)=1;%ODEOPTIONS(11:12)=[-1,-1];
 
405
 
 
406
deff('jj=j1(t,x)','jj=A')
 
407
 
 
408
xnotband=ode('st',x0,t0,t,f,j1);
 
409
 
 
410
 
 
411
 
 
412
%ODEOPTIONS(6)=4;//banded jacobian external given
 
413
 
 
414
%ODEOPTIONS(11:12)=[3,2];
 
415
 
 
416
deff('jj=j(t,x)','jj=J')
 
417
 
 
418
xband=ode('st',x0,t0,t,f,j);
 
419
 
 
420
if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 
421
 
 
422
 
 
423
%ODEOPTIONS(6)=5;//banded jacobian evaluated
 
424
 
 
425
%ODEOPTIONS(11:12)=[3,2];
 
426
 
 
427
deff('jj=j(t,x)','jj=J')
 
428
 
 
429
xband=ode('st',x0,t0,t,f,j);
 
430
 Warning: Jacobian external is given, but 
 
431
 not used!,  see %ODEOPTIONS(6)
 
432
 
 
433
 
 
434
if norm(xnotband-xband,1)>Leps then bugmes();quit;end
 
435
 
 
436
 
 
437
 
 
438
//            Test of %ODEOPTIONS(7)
 
439
 
 
440
//%ODEOPTIONS(7)=mxstep  ---> maximum number od steps allowed
 
441
 
 
442
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
443
 
 
444
//provisional values as default
 
445
 
 
446
jacflag=2;ml=-1;mu=-1;
 
447
 
 
448
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
449
 
 
450
%ODEOPTIONS(7)=10;
 
451
 
 
452
//ode(x0,t0,t,f);  // ---> Non convergence
 
453
 
 
454
 
 
455
//            Test of %ODEOPTIONS(8:9)
 
456
 
 
457
//%ODEOPTIONS(8:9)=[maxordn,maxords] ---> maximum order for nonstiff and stiff
 
458
 
 
459
itask=1;tcrit=0;h0=0;hmax=0;hmin=0;ixpr=0;mxstep=0;maxordn=0;maxords=0;
 
460
 
 
461
//provisional values as default
 
462
 
 
463
jacflag=2;ml=-1;mu=-1;
 
464
 
 
465
%ODEOPTIONS=[itask,tcrit,h0,hmax,hmin,jacflag,mxstep,maxordn,maxords,ixpr,ml,mu];
 
466
 
 
467
%ODEOPTIONS(8:9)=[0,0]; //---> default values
 
468
 
 
469
wref=ode(x0,t0,t,f); //just for computing reference
 
470
 
 
471
 
 
472
%ODEOPTIONS(8:9)=[4,3];
 
473
 
 
474
ww=ode(x0,t0,t,f);norm(wref-ww,1);
 
475
 
 
476
 
 
477
%ODEOPTIONS(8:9)=[12,5];
 
478
 
 
479
if norm(wref-ode(x0,t0,t,f),1)>Leps then bugmes();quit;end
 
480
 
 
481
 
 
482
//using stiff method
 
483
 
 
484
 
 
485
%ODEOPTIONS(9)=0;
 
486
 
 
487
wref=ode('st',x0,t0,t,f);
 
488
 
 
489
%ODEOPTIONS(9)=5;
 
490
 
 
491
if norm(wref-ode('st',x0,t0,t,f),1) >Leps then bugmes();quit;end //=0
 
492
 
 
493
%ODEOPTIONS(9)=4;
 
494
 
 
495
if norm(wref-ode('st',x0,t0,t,f),1)  >Leps then bugmes();quit;end  //small
 
496
 
 
497
 
 
498
 
 
499
//using nonstiff method
 
500
 
 
501
 
 
502
%ODEOPTIONS(8)=0;
 
503
 
 
504
wref=ode('ad',x0,t0,t,f);
 
505
 
 
506
%ODEOPTIONS(8)=12;
 
507
 
 
508
if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end  //=0
 
509
 
 
510
%ODEOPTIONS(8)=5;
 
511
 
 
512
if norm(wref-ode('ad',x0,t0,t,f),1) >Leps then bugmes();quit;end   //small
 
513
 
 
514
 
 
515
//mixed
 
516
 
 
517
%ODEOPTIONS(8:9)=[5,12];
 
518
 
 
519
wref=ode(x0,t0,t,f);
 
520
 
 
521
%ODEOPTIONS(8:9)=[4,10];
 
522
 
 
523
if norm(ode(x0,t0,t,f)-wref,1)>Leps then bugmes();quit;end   //small
 
524
 
 
525
 
 
526
 
 
527
A=diag([-10,-0.01,-1]);
 
528
 
 
529
deff('uu=u(t)','uu=sin(t)');
 
530
 
 
531
B=rand(3,1);
 
532
 
 
533
deff('y=f(t,x)','y=A*x+B*u(t)')
 
534
Warning :redefining function: f                       
 
535
 
 
536
 
 
537
%ODEOPTIONS(1)=2;
 
538
 
 
539
yy1=ode('stiff',[1;1;1],0,1,f);
 
540
 
 
541
yy2=ode('stiff',[1;1;1],0,2,f);
 
542
 
 
543
%ODEOPTIONS(1)=3;
 
544
 
 
545
yy1=ode('stiff',[1;1;1],0,1,f);
 
546
 
 
547
yy2=ode('stiff',[1;1;1],0,2,f);
 
548
 
 
549
 
 
550
clear %ODEOPTIONS;
 
551
 
 
552
rand('seed',0);rand('normal');
 
553
 
 
554
nx=20;A=rand(nx,nx);A=A-4.5*eye();
 
555
 
 
556
deff('y=f(t,x)','y=A*x')
 
557
Warning :redefining function: f                       
 
558
 
 
559
 
 
560
deff('J=j(t,x)','J=A')
 
561
Warning :redefining function: j                       
 
562
 
 
563
 
 
564
//%ODEOPTIONS(1)=1;
 
565
 
 
566
y2=ode('stiff',ones(nx,1),0,2,f,j);
 
567
 
 
568
[y1,w,iw]=ode('stiff',ones(nx,1),0,1,f,j);
 
569
 
 
570
y2p=ode('stiff',y1,1,2,f,j,w,iw);
 
571
 
 
572
y12=ode('stiff',ones(nx,1),0,[1,2],f,j);
 
573
 
 
574
norm(y12(:,2)-y2p);
 
575
 
 
576
yaf=ode('adams',ones(nx,1),0,2,f,j);
 
577
 
 
578
yaj=ode('adams',ones(nx,1),0,2,f,j);
 
579
 
 
580
ysf=ode('stiff',ones(nx,1),0,2,f,j);
 
581
 
 
582
ysj=ode('stiff',ones(nx,1),0,2,f,j);
 
583
 
 
584
 
 
585
deff('xd=f(t,x)','xd=A*x+B*sin(3*t)')
 
586
Warning :redefining function: f                       
 
587
 
 
588
 
 
589
A=rand(10,10)-4.5*eye();B=rand(10,1);
 
590
 
 
591
x=ode(ones(10,1),0,[1,2,3],f);
 
592
 
 
593
//link('fexab.o','fexab')
 
594
 
 
595
if norm(x-ode(ones(10,1),0,[1,2,3],'fexab'),1) > 1.d-10 then bugmes();quit;end
 
596
 
 
597