StepperRoss Example Problem


MPD78
08-12-2009, 07:49 AM
Hello all,

In the example problem that uses stepperross.h, on page 940 of the online book, I am slightly confused by the struct rhs.

My confusion is this;

What happend to the following two lines?

Doub eps;
rhs(Doub epss) : eps(epss) {}

If I code the example verbatim I recieve an error upon compilation that states

"error C2440: 'initializing' : cannot convert from 'double' to 'rhs'

Now if I add the two (assumed) missing lines into the struct rhs, and set nvar = 3, the program compilies correct
and the program executes correctly with 11 good steps and 0 bad steps

Here is my main() function

#include "odeint.h"
#include "stepperdopr5.h"
#include "stepperbs.h
#include "stepperross.h"

struct rhs {
// Code from page 940 of the online text
}

void jacobian .... {
// code from page 940 of the online text

}


int main()
{
// Code to open up output file

const Int nvar = 3;
const Doub atol=1.0e-5,rtol=atol, h1=2.9e-4, hmin=0.0, x1=0.0, x2=50.0;
VecDoub ystart(nvar);
ystart[0]=1;
ystart[1]=1;
ystart[2]=0;
Output out(-1);
rhs d(1.0e-5);
Odeint<StepperRoss<rhs> > ode(ystart,x1,x2,atol,rtol,h1,hmin,out,d);
ode.integrate();

// Code to write data to output file

return 0;
}



Results

Example Results

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

0 1 1 0
0.00029 1 1 -2.367e-006
0.00203 0.99998 1 -3.613e-006
0.0078883 0.99993 1.0001 -3.6974e-006
0.042023 0.99961 1.0004 -3.7112e-006
0.24683 0.99771 1.0023 -3.7021e-006
1.4757 0.98634 1.0137 -3.6422e-006
8.8488 0.9194 1.0806 -3.3009e-006
21.909 0.80713 1.1929 -2.7691e-006
35.318 0.70144 1.2986 -2.3098e-006
49.71 0.59958 1.4004 -1.9008e-006
50 0.59766 1.4023 -1.8934e-006

Number of good steps 11 The number of bad steps 0
So my overall question is;

Should the following two lines be be added to the struct rhs?

Doub eps;
rhs(Doub epss) : eps(epss) {}

Or am I way off in the wrong direction?

Thanks
Matt

Idiot Disclaimer

"I am an engineer not a highly trained programmer. So please be gentle."

Saul Teukolsky
08-12-2009, 10:14 AM
Hi Matt,

Unlike the example on p.905-906, this example has no parameter like eps in the right-hand side. So you don't need the two lines you are concerned about. Instead of the line
rhs d(1.0e-5);
in the calling program (which does nothing since 1.0e-5 is not used anywhere), you just need
rhs d;
This invokes the default constructor in rhs (which you don't need to supply), and everything will work.

Best,
Saul Teukolsky

MPD78
08-12-2009, 10:19 AM
Saul,

Thanks.

I changed the code accordingly and all works fine with identical results as those posted above.

Matt