Hello folks,

so I am still trying to work with the NR. Now I want to solve ODEs via chapter 17.
I managed to get the quick-example startet (§ 17.0.4) but I have no idead what the graph shall tell me.
So I want to solve some trivial ODEs where I know solution like (radioactive) decay or the harmonic oscillator.

But again I have no idea how to use the NR routines. In the book there are some explanations about the idea of solving an ODEs and how an algorithmen works in principle but I am missing some connection.

So let's take the decay, which is simply \cdot x = - c * x where c is simply a constant and x the number of anything decaying e.g. radioactive material.

Now I have to rewrite the formula \Delta x = -c * x * \Delta t.

And now? NR says I have to include it in a constructor (class):

class decay
{
public:
};


But now I have no idea how to proceed! :-(
I cannot simply adapt the quick-example for this case, since it is another sort of ODE.

Sincerely,
Quantum.

Oh,

it seems, that this is the wrong Thread?

Maybe an admin can shift it to Methods: Chapters 16 and 17.

As an user I am not allowed to erase my own post (and place it elsewhere).

Cheers,
Quantum.

According to the owners, for now, all Version 3 questions are in this part of the forum. The Methods: Chapters 16 and 17 part is for previous (obsolete) versions.

So, dis must be de place.



\dot x = - c x

Odeint is used to solve a system of first-order equations. The right-hand side is calculated with a functor (class with an overloaded () operator)

If you have a single first order equation, you use vectors with size equal to 1.


#include "nr3.h"
#include "stepper.h"
#include "stepperdopr5.h"
#include "odeint.h"

//
// Solve a single first-order ode
// dy/dx = -c*x
//
// This class creates a functor to be used for this first order ODE
//

struct rhs_van
{
Doub c;
Doub eps;

//
// Constructor for a particular value of c
//
// The "eps" argument is optional. If none is
// given, the default value will be used.
//
rhs_van(Doub cc, Doub epss = 1.0e-5) : c(cc), eps(epss) {}

//
// This returns the value of dy/dx for given x and y.
// Note that this particular dydx doesn't depend on x
// but this is the general form required by odeint
//
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0]= -c * y[0];
}
};

//
// We picked an easy one for which the answer is known
//
Doub actual(Doub y0, Doub c, Doub x)
{
return y0 * exp(-c * x);
}

int main()
{
// This is the constant used in dy/dx
Doub c = 1.0;

Doub eps = 1.0e-4;
Doub h1 = 0.1;
Doub hmin = 0.0;

Doub atol = 1.0e-5; // absolute tolerance
Doub rtol = atol/2.0; // relative tolerance

Doub x0 = 0.0; // Start of interval
Doub x1 = 1.0; // End of interval
Doub y0 = 1.0; // Value of dy/dx at x0

// A vector used for initial value
VecDoub ystart(1);
ystart[0] = y0;

// Output will be calculated for 20 points after x0
Output out(20);

// Here's what odeint will use to calculate the right hand side
rhs_van d(c, eps);

// Finally: here's the Odeint object
Odeint<StepperDopr5<rhs_van> >ode(ystart, x0, x1, atol, rtol, h1, hmin, out, d);

// Do the work of calculating output values
ode.integrate();

//
// Print values at each output point
//
cout << scientific << setprecision(5);
cout << " x Approx Actual" << endl;
cout << " --------------------------------------" << endl;
for (Int i = 0; i < out.count; i++) {
cout << setw(13) << out.xsave[i]
<< setw(13) << out.ysave[0][i]
<< setw(13) << actual(y0, c, out.xsave[i])
<< endl;
}
return 0;
}


Output:

x Approx Actual
--------------------------------------
0.00000e+00 1.00000e+00 1.00000e+00
5.00000e-02 9.51229e-01 9.51229e-01
1.00000e-01 9.04837e-01 9.04837e-01
1.50000e-01 8.60707e-01 8.60708e-01
2.00000e-01 8.18728e-01 8.18731e-01
2.50000e-01 7.78796e-01 7.78801e-01
3.00000e-01 7.40814e-01 7.40818e-01
3.50000e-01 7.04685e-01 7.04688e-01
4.00000e-01 6.70319e-01 6.70320e-01
4.50000e-01 6.37629e-01 6.37628e-01
5.00000e-01 6.06532e-01 6.06531e-01
5.50000e-01 5.76950e-01 5.76950e-01
6.00000e-01 5.48811e-01 5.48812e-01
6.50000e-01 5.22044e-01 5.22046e-01
7.00000e-01 4.96584e-01 4.96585e-01
7.50000e-01 4.72366e-01 4.72367e-01
8.00000e-01 4.49329e-01 4.49329e-01
8.50000e-01 4.27416e-01 4.27415e-01
9.00000e-01 4.06571e-01 4.06570e-01
9.50000e-01 3.86743e-01 3.86741e-01
1.00000e+00 3.67881e-01 3.67879e-01


Regards,

Dave

Uh,
thank you for the explanation. This is exactly what I am looking for and missing in the NR.
Now, according to "every answer leads to ten new questions" I have some questions. :)

1.
Ho do you manage to write the DGL? Is there some sort of LaTex Plug-In?

2.
What is eps?
First I looked at the Van-der Pol ODE and there it is a parameter for the oscillating. But for a simple decay ODE I do not need an oscillating-parameter! :confused:

3.
Where do I get an explanation for atol and rtol? They are introduced in §17.2 with some remarks but what do they mean? Are they important for ODEs or at all for numerical calculations?

4.
For the initial values you write:

Doub y0 = 1.0; // Value of dy/dx at x0

I think the comment should be "Value of y0 at x0"?

Thanks in advance.

Best regards,
Quantum.

1.
Ho do you manage to write the DGL?It's magic. See Footnote.

Here's the drill:
When composing a post, put some TeX stuff on a line by itself. Put at the beginning of the line and put at the end of the line. Use the "Preview Post" button to see if it looks like you want it to look.

2.
What is eps?
Sorry. That was a sloppy holdover from a function I had tested that had an eps parameter. Since the function in your problem does not have an eps anywhere in the formula, that parameter is completely superfluous.

So the constructor looks like:

//
// Constructor for a particular value of c
//
rhs_van(Doub cc) : c(cc) {}


And here's how you create the object:

// Here's what odeint will use to calculate the right hand side
rhs_van d(c);



3.
Where do I get an explanation for atol and rtol?It's an adaptive algorithm that adjusts step sizes to get to each calculated output point. The values of atol and rtol give the maximum allowable (estimated) absolute error and maximum allowable (estimated) relative error at each point. If you make them too big, the truncation error may be too large. If you make them unreasonably small, it may require too many steps to be practical (depending on the function, obviously).


4.
For the initial values you write:

Doub y0 = 1.0; // Value of dy/dx at x0

Indeed, the comment should be "Value of y for x = x0" or some such thing.


Regards,

Dave

Footnote:
"Any sufficiently advanced technology
is indistinguishable from magic."
---Arthur C. Clarke

Thank you very much! :)

I think more questions will follow. :D

Have a nice weekend! :)

Best regards,
Quantum.

So I am back.
As a physician I am first of all interested in solving the harmonic oscillator equation:

\ddot x + \omega^{2}x = 0


Now I found two answers for my question "how to" ?

First, what NR says for itself: Reduce the equation to a set of two 1st order equations.
I found example (http://www.kramann.info/85_Archiv_Simulationstechnik/02_Newton/01_LinearSchwinger/index.php) (sorry, in german!).
But now, how to coupe with that matrix?

Second,
NR says explicitly that you can use Stoermer's rule (§17.4) for second order ODE if there is no first derivate on the right side. That would fit with my first trivial for oscillating without any dumping or external force.

Thank you very much for your patience!

Best regards,
Quantum.

...
I found example (http://www.kramann.info/85_Archiv_Simulationstechnik/02_Newton/01_LinearSchwinger/index.php)

//
// Demo of StepperDopr5 with Odeint to solve 2nd order ODE
//
// davekw7x
//
#include "nr3.h"
#include "stepper.h"
#include "stepperdopr5.h"
#include "odeint.h"

// Some versions of cmath define M_PI, some do not
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

int numcalc = 0;

//
// Second order ode: x dotdot = -C/m x
//
// Let v = xdot so vdot = x dotdot = -C/m x
//
// Then we solve equivalent system of first-order odes
// xdot = v (We want to solve this for x)
// vdot = -C/m * x
//
// Fit these two into a vector y for solution by Ode:
// y[0] dot = y[1] (We want to solve this for y[0])
// y[1] dot = -C/m * y[0]
//
// Here's a functor to be used for the system of first order ODEs


struct rhs_van
{
Doub c; // Spring constant
Doub m; // Mass

//
// Constructor for a particular value of c
//
// Note that this particular dydx doesn't depend on x
// but this is the general form required by odeint
//
rhs_van(Doub cc, Doub mm) : c(cc),m(mm) {}

//
// This returns the values of dy/dx
//
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0] = y[1]; // y dot
dydx[1] = -(c/m) * y[0]; // y dotdot
++numcalc;
}
};

//
Doub actual(Doub y0, Doub y0p, Doub c, Doub m, Doub x)
{
Doub omega = sqrt(c/m);
return y0 * cos(omega*x) + y0p/omega * sin(omega*x);
}

int main()
{
// This is the constant used in dy/dx
Doub c = 4.0;
Doub m = 1.0;

Doub h1 = 0.001;
Doub hmin = 0.0;

Doub atol = 1.0e-7; // absolute tolerance
Doub rtol = atol/2.0; // relative tolerance

Doub x0 = 0.0; // Start of interval
Doub x1 = M_PI; // End of interval
Doub y0 = 1.0; // Value of y at x0
Doub y0p = 1.0; // Value of dy/dx at x0

// A vector used for initial value
VecDoub ystart(2);
ystart[0] = y0;
ystart[1] = y0p;

// Output will be calculated for 20 points after x0
Output out(20);

// Here's what odeint will use to calculate the right hand side
rhs_van d(c, m);

// Finally: here's the ode object
Odeint<StepperDopr5<rhs_van> > ode(ystart, x0, x1, atol, rtol, h1, hmin, out, d);

// Do the work of calculating output
ode.integrate();

//
// Print values at each output point
//
cout << "c = " << c << ", m = " << m << ", omega = " << sqrt(c/m)
<< endl << endl;
cout << scientific << setprecision(5);
cout << " x Approx Actual" << endl;
cout << " --------------------------------------" << endl;
for (int i = 0; i < out.count; i++) {
cout << setw(13) << out.xsave
<< setw(13) << out.ysave[0][i]
<< setw(13) << actual(y0, y0p, c, m, out.xsave[i])
<< endl;
}

cout << endl << "Number of calculations = " << numcalc << endl;
return 0;
}


Output:

c = 4, m = 1, omega = 2

x Approx Actual
--------------------------------------
0.00000e+00 1.00000e+00 1.00000e+00
1.57080e-01 1.10557e+00 1.10557e+00
3.14159e-01 1.10291e+00 1.10291e+00
4.71239e-01 9.92294e-01 9.92294e-01
6.28319e-01 7.84545e-01 7.84545e-01
7.85398e-01 5.00000e-01 5.00000e-01
9.42478e-01 1.66511e-01 1.66511e-01
1.09956e+00 -1.83277e-01 -1.83277e-01
1.25664e+00 -5.15124e-01 -5.15124e-01
1.41372e+00 -7.96548e-01 -7.96548e-01
1.57080e+00 -1.00000e+00 -1.00000e+00
1.72788e+00 -1.10556e+00 -1.10557e+00
1.88496e+00 -1.10291e+00 -1.10291e+00
2.04204e+00 -9.92294e-01 -9.92294e-01
2.19911e+00 -7.84545e-01 -7.84545e-01
2.35619e+00 -5.00000e-01 -5.00000e-01
2.51327e+00 -1.66511e-01 -1.66511e-01
2.67035e+00 1.83277e-01 1.83277e-01
2.82743e+00 5.15124e-01 5.15124e-01
2.98451e+00 7.96548e-01 7.96548e-01
3.14159e+00 1.00000e+00 1.00000e+00

Number of calculations = 259


Second,
...you can use [i]Stoermer's rule (§17.4) for second order ODE if there is no first derivat[ive] on the right side....

//
// Demo of StepperStoerm with Odeint to solve 2nd order conservative ODE
//
// davekw7x
//
#include "nr3.h"
#include "stepper.h"
#include "stepperbs.h"
#include "stepperstoerm.h"
#include "odeint.h"

// Some versions of cmath define M_PI, some do not
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

int numcalc = 0;

//
// Second order ode: x dotdot = -C/m x
//

struct rhs_van
{
Doub c; // Spring constant
Doub m; // Mass

//
// Constructor for a particular value of c and m
//
rhs_van(Doub cc, Doub mm) : c(cc),m(mm) {}

//
// This returns the value of d^2y/dx^2 for given x and y
// Note that for this ode, dy^2/dx^2 doesn't depend on x
// but this is the general form required by odeint
//
//
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0] = -(c/m) * y[0];
++numcalc;
}
};

Doub actual(Doub y0, Doub y0p, Doub c, Doub m, Doub x)
{
Doub omega = sqrt(c/m);
return y0 * cos(omega*x) + y0p/omega * sin(omega*x);
}


int main()
{
Doub c = 4.0; // Spring constant
Doub m = 1.0; // Mass

Doub h1 = 0.001;
Doub hmin = 0.0;

Doub atol = 1.0e-7; // absolute tolerance
Doub rtol = atol/2.0; // relative tolerance

Doub x0 = 0.0; // Start of interval
Doub x1 = M_PI; // End of interval
Doub y0 = 1.0; // Value of y at x0
Doub y0p = 1.0; // Value of dy/dx at x0

// A vector used for initial value
VecDoub ystart(2);
ystart[0] = y0;
ystart[1] = y0p;

// Output will be calculated for 20 points after x0
Output out(20);

// Here's what odeint will use to calculate the right hand side
rhs_van d(c, m);

// Finally: here's the ode object
Odeint<StepperStoerm<rhs_van> > ode(ystart, x0, x1, atol, rtol, h1, hmin, out, d);

// Do the work of calculating output
ode.integrate();

//
// Print values at each output point
//
cout << "c = " << c << ", m = " << m << ", omega = " << sqrt(c/m)
<< endl << endl;
cout << scientific << setprecision(5);
cout << " x Approx Actual" << endl;
cout << " --------------------------------------" << endl;
for (int i = 0; i < out.count; i++) {
cout << setw(13) << out.xsave[i]
<< setw(13) << out.ysave[0][i]
<< setw(13) << actual(y0, y0p, c, m, out.xsave[i])
<< endl;
}

cout << endl << "Number of calculations = " << numcalc << endl;
return 0;
}


Output


c = 4, m = 1, omega = 2

x Approx Actual
--------------------------------------
0.00000e+00 1.00000e+00 1.00000e+00
1.57080e-01 1.10557e+00 1.10557e+00
3.14159e-01 1.10291e+00 1.10291e+00
4.71239e-01 9.92294e-01 9.92294e-01
6.28319e-01 7.84545e-01 7.84545e-01
7.85398e-01 5.00000e-01 5.00000e-01
9.42478e-01 1.66511e-01 1.66511e-01
1.09956e+00 -1.83277e-01 -1.83277e-01
1.25664e+00 -5.15124e-01 -5.15124e-01
1.41372e+00 -7.96548e-01 -7.96548e-01
1.57080e+00 -1.00000e+00 -1.00000e+00
1.72788e+00 -1.10557e+00 -1.10557e+00
1.88496e+00 -1.10291e+00 -1.10291e+00
2.04204e+00 -9.92294e-01 -9.92294e-01
2.19911e+00 -7.84545e-01 -7.84545e-01
2.35619e+00 -5.00000e-01 -5.00000e-01
2.51327e+00 -1.66511e-01 -1.66511e-01
2.67035e+00 1.83277e-01 1.83277e-01
2.82743e+00 5.15124e-01 5.15124e-01
2.98451e+00 7.96548e-01 7.96548e-01
3.14159e+00 1.00000e+00 1.00000e+00

Number of calculations = 148



Regards,

Dave

Thank you again a lot Dave.

What I wonder now:
How I will proceed in the future? I mean, I cannot always ask for an example.

Are my C++ programming skills just insufficient so far to use NR properly? So I will just get used to it?
Or is there something like a practical guide to NR?

And again some questions! :)

First, one of the differences between the decay ode \dot x = -c x and the harmonic oscillator \ddot x + \omega^{2} x = 0 is the the line about A vector used for initial value


// A vector used for initial value
VecDoub ystart(2);
ystart[0] = y0;
ystart[1] = y0p;


According to the owners, for now, all Version 3
// A vector used for initial value
VecDoub ystart(1);
ystart[0] = y0;


I wonder about the number in the term ystart(NUMBER). Is this number important to give the number of boundary conditions?

Second if you define the operator for the harmonic oscillator (stepper method) I do not understand what the ++numcalc is for?

void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0] = y[1]; // y dot
dydx[1] = -(c/m) * y[0]; // y dotdot
++numcalc;
}


Best regards!

Quantum.

...is there something like a practical guide to NRMaybe you can find a textbook or an on-line tutorial on C++ that you can work with. Some people like this one: http://www.cplusplus.com/doc/tutorial/. Then look at NR source code. The listings in the text have additional helpful comments but it's not always easy to read and understand Other People's Code, even with good documentation. Some people have a knack for it; for most of us it's a learned skill.

the number in the term ystart(NUMBER)

For StepperDopr5:
The vector is declared with a size equal to the number of equations. If the original ODE is of order n, there will be n first order ODEs and the size of the initial value vector will be n. Then it is populated with initial conditions using array-style indexed notation (the [] stuff).


For StepperStoerm:
The vector is declared with a size equal to the order of the equation. If the order is n, there will be n initial conditions (the value of the variable and the values of the first n-1 derivatives). Therefore, the vector is declared with size equal to n.

...++numcalc...That's just so that the program can report how many times the functor was called, and has nothing whatsoever to do with the solution itself. Sometimes people are interested in comparing methods to see whether one is more "efficient" than another or some such thing.


Regards,

Dave

Thanks.
This explains a lot.

So for the next step I thought it would be necessary to solve partial differential equations. But after reading some pages I have the impression that this is a big step from ode to pde?
The NR say that it is possible to fill whole volumes about solving pdes.

I googled it and found out: Solving PDEs with C++ (http://www.amazon.com/Solving-PDEs-Computational-Science-Engineering/dp/0898716012)

This seems to be some sort of standard?

Sorry for beeing a little off topic.

Best regards,
Quantum.

Edit: So now I am on vacation until 29th august. I will respond after my return.

According to the owners, for now, all Version 3 questions are in this part of the forum. The Methods: Chapters 16 and 17 part is for previous (obsolete) versions.

So, dis must be de place.




Odeint is used to solve a system of first-order equations. The right-hand side is calculated with a functor (class with an overloaded () operator)

If you have a single first order equation, you use vectors with size equal to 1.


#include "nr3.h"
#include "stepper.h"
#include "stepperdopr5.h"
#include "odeint.h"

//
// Solve a single first-order ode
// dy/dx = -c*x
//
// This class creates a functor to be used for this first order ODE
//

struct rhs_van
{
Doub c;
Doub eps;

//
// Constructor for a particular value of c
//
// The "eps" argument is optional. If none is
// given, the default value will be used.
//
rhs_van(Doub cc, Doub epss = 1.0e-5) : c(cc), eps(epss) {}


// This returns the value of dy/dx for given x and y.
// Note that this particular dydx doesn't depend on x
// but this is the general form required by odeint
//
void operator() (const Doub x, VecDoub_I &y, VecDoub_O &dydx)
{
dydx[0]= -c * y[0];
}
};

//
// We picked an easy one for which the answer is known
//
Doub actual(Doub y0, Doub c, Doub x)
{
return y0 * exp(-c * x);
} // التسويق (http://www.nile7.com/readings/التسويق/)

int main()
{
// This is the constant used in dy/dx
Doub c = 1.0;

Doub eps = 1.0e-4;
Doub h1 = 0.1;
Doub hmin = 0.0;

Doub atol = 1.0e-5; // absolute tolerance
Doub rtol = atol/2.0; // relative tolerance

Doub x0 = 0.0; // Start of interval
Doub x1 = 1.0; // End of interval
Doub y0 = 1.0; // Value of dy/dx at x0

// A vector used for initial value
VecDoub ystart(1);
ystart[0] = y0;
اشهار منتديات (http://www.nile7.com)
// Output will be calculated for 20 points after x0
Output out(20);

// Here's what odeint will use to calculate the right hand side
rhs_van d(c, eps);

// Finally: here's the Odeint object
Odeint<StepperDopr5<rhs_van> >ode(ystart, x0, x1, atol, rtol, h1, hmin, out, d);

// Do the work of calculating output values
ode.integrate();

//
// Print values at each output point
//
cout << scientific << setprecision(5);
cout << " x Approx Actual" << endl;
cout << " --------------------------------------" << endl;
for (Int i = 0; i < out.count; i++) {
cout << setw(13) << out.xsave[i]
<< setw(13) << out.ysave[0][i]
<< setw(13) << actual(y0, c, out.xsave[i])
<< endl;
}
return 0;
}


Output:

x Approx Actual
--------------------------------------
0.00000e+00 1.00000e+00 1.00000e+00
5.00000e-02 9.51229e-01 9.51229e-01
1.00000e-01 9.04837e-01 9.04837e-01
1.50000e-01 8.60707e-01 8.60708e-01
2.00000e-01 8.18728e-01 8.18731e-01
2.50000e-01 7.78796e-01 7.78801e-01
3.00000e-01 7.40814e-01 7.40818e-01
3.50000e-01 7.04685e-01 7.04688e-01
4.00000e-01 6.70319e-01 6.70320e-01
4.50000e-01 6.37629e-01 6.37628e-01
5.00000e-01 6.06532e-01 6.06531e-01
5.50000e-01 5.76950e-01 5.76950e-01
6.00000e-01 5.48811e-01 5.48812e-01
6.50000e-01 5.22044e-01 5.22046e-01
7.00000e-01 4.96584e-01 4.96585e-01
7.50000e-01 4.72366e-01 4.72367e-01
8.00000e-01 4.49329e-01 4.49329e-01

8.50000e-01 4.27416e-01 4.27415e-01
9.00000e-01 4.06571e-01 4.06570e-01
9.50000e-01 3.86743e-01 3.86741e-01
1.00000e+00 3.67881e-01 3.67879e-01


Regards,

Dave

very good work