This is what I think is an error in StepperSie's implementation (Webnote 24, rev 1).

In function 'step', there are several static variables, some of which are initialized there (e g, first_step). When the integrator is called the first time, these variables are initialized to their proper values. But if the integrator is called a second time, the variables are not initialized, and so do not have appropriate values (e g, first_step will be false, when it should be true).

The fix for this seems to be having an 'initializeMe' argument to step, which causes 'step' to do these initializations, or a similar value in the base class, accessible to the applications that are using the ODE solver. Alternatively, these variables could be moved into that class, so the user could set them directly, although this puts more burden on the user. Best to have everything happen automatically.

Hi Woody,

I'm not sure there is an error. If you call the integrator twice, with two separate invocations of the constructor like the following, there's no problem:
Odeint<StepperSie<rhs> > integ(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
integ.integrate();
...
//choose new values of x1, x2, etc
Odeint<StepperSie<rhs> > integ2(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
integ2.integrate();

I agree that if you just have one constructor, and then try
integ.integrate();
...
integ.integrate();
where you simply reset some of the variables between calls, you are likely to have trouble. Can you give an example of what you would like to do that you believe doesn't currently work?

Best,
Saul Teukolsky

Hello all,

I tried the method you described above and this is the output I recieved. I used the stiff ODE example given in the book. First I integrated the system from x=0 to x=50 and then I integrated the system from x=0 to x=100. There is something odd about the first time step after x=50 in the second integration. The time step becomes zero again for some reason.

Results for x = 0 to x = 50

Integration from x=0 to x=50
============================================
xsave[i] ysave[0][i] ysave[1][i]

0 1 1
0.00029 0.999997 1
0.00230753 0.999978 1.00002
0.0163435 0.999847 1.00015
0.113992 0.998941 1.00106
0.793335 0.992644 1.00735
3.9797 0.963349 1.03665
10.8445 0.901699 1.0983
25.6343 0.776746 1.22325
37.8172 0.682876 1.31712
50 0.597656 1.40234

Results for x = 0 to x = 100

Integration from x=0 to x=100
============================================

xsave[i] ysave[0][i] ysave[1][i]

0.000000 1.000000 1.000000
0.000290 0.999997 1.000001
0.002308 0.999978 1.000019
0.016344 0.999847 1.000149
0.113992 0.998941 1.001056
0.793335 0.992644 1.007352
3.979699 0.963349 1.036647
10.844513 0.901699 1.098297
25.634306 0.776746 1.223251
37.817153 0.682876 1.317122
50.000000 0.597656 1.402342
0.000000 0.597656 1.402342
0.000290 0.597654 1.402344
0.002308 0.597641 1.402357
0.016344 0.597548 1.402450
0.113992 0.596900 1.403098
0.793335 0.592408 1.407590
3.246984 0.576410 1.423588
11.176324 0.527086 1.472913
22.504326 0.462756 1.537243
32.890024 0.409824 1.590174
40.464289 0.374658 1.625341
48.038554 0.342220 1.657779
55.391773 0.313189 1.686810
62.744991 0.286434 1.713565
70.165853 0.261594 1.738405
77.586714 0.238779 1.761220
85.173762 0.217399 1.782600
92.760809 0.197842 1.802157
100.000000 0.180754 1.819246

Why does the timestep go to zero after 50?

With one call to StepperSie the results from integrating from x=0 to x=100 are not the same when using two calls to StepperSie.

Integration from x=0 to x=100
============================================
xsave[i] ysave[0][i] ysave[1][i]

0 1 1
0.00029 0.999997 1
0.00230753 0.999978 1.00002
0.0163435 0.999847 1.00015
0.113992 0.998941 1.00106
0.793335 0.992644 1.00735
3.9797 0.963349 1.03665
10.8445 0.901699 1.0983
25.6343 0.776746 1.22325
28.1533 0.75664 1.24336
30.6724 0.736893 1.2631
36.0995 0.69559 1.30441
47.7919 0.612455 1.38754
58.188 0.54525 1.45475
68.5841 0.484218 1.51578
79.9276 0.424352 1.57565
91.271 0.371075 1.62892
100 0.334244 1.66576

Can anyone explain this difference?

Thanks
Matt

Hi Matt,

You need to create a new Output object as well, otherwise you are just adding the new output values to the old ones (your results for 0 to 100 are just the results for 0 to 50 followed by the actual results for 0 to 100). So before the second integration you could declare
Output out2(20);
and pass out2 as an argument to integ2. (Or you could use the original out object but print only the new saved steps.)

Saul Teukolsky