Hello,

When buying the source code from Numerical Receipes (NR), some example code was also included. The example code "xvoltra.cpp" is an example on how to solve the Volterra equation. This file is attached here ( I had to rename it to "xvoltra.txt" ).

In this example code given by NR we have N=50, where N = number of steps. This works just fine. Now, running for larger number of steps, e.g. N=500, the numerical solution divergates. Even for very small stepsize "H" this is still a problem. See attached PDF-file.

Why do I always run into divergence problems, even when applying NR's own example code?

Regards, "perr"

Are you referencing the NR3 book ? There is no example file xvoltra.cpp in that book.

There is an example in section 19.3.2.

Basically the more descriptive you can be with your question the better the chance you have of recieving help on this forum.

Thanks
Matt

Thank you for your reply.
Now I have updated my question. This time I am more descriptive. In addition, I have attached the file "xvoltra.txt".

Can anybody help me?: Why do I always run into divergence problems, even when applying NR's own example code?

Best, "perr"

Your attached text file is in very hard to read format. Can you just post the code inside of a code block?

It is tool icon that looks like the # sign in the thread tools.

Thanks
Matt

Hi again,

I don't understand the issue concerning read format. I therefore put the code for "xvoltra.cpp" here:

#include <iostream>
#include <iomanip>
#include <cmath>
#include "nr.h"
using namespace std;

// Driver for routine voltra

DP g(const int k, const DP t)
{
return (k == 0 ? cosh(t)+t*sin(t) :
2.0*sin(t)+t*(SQR(sin(t))+exp(t)));
}

DP ak(const int k, const int l, const DP t, const DP s)
{
return ((k == 0) ?
(l == 0 ? -exp(t-s) : -cos(t-s)) :
(l == 0 ? -exp(t+s) : -t*cos(s)));
}

int main(void)
{
const int N=50,M=2;
const DP H=0.01;
int nn;
DP t0=0.0;
Vec_DP t(N);
Mat_DP f(M,N);

NR::voltra(t0,H,t,f,g,ak);
// exact soln is f[1]=exp(-t), f[2]=2sin(t)
cout << " abscissa voltra real";
cout << " voltra real" << endl;
cout << " answer1 answer1";
cout << " answer2 answer2" << endl << endl;
cout << fixed << setprecision(6);
for (nn=0;nn<N;nn++) {
cout << setw(12) << t[nn] << setw(13) << f[0][nn];
cout << setw(13) << exp(-t[nn]) << setw(13) << f[1][nn];
cout << setw(13) << 2.0*sin(t[nn]) << endl;
}
return 0;
}

This is the code given by NR themselves.

Why do I always run into divergence problems, even when applying NR's own example code?

Best, "perr"

Okay, thanks for posting the code.

I am not the best with these functions but I have recieved some good help here and I am trying to return it. Maybe you can help me now.

Below is my version of your code. See the comments inside the code.

My code isn't working properly. (This is where you can help me.) I keep getting an error that says the function g does not evaluate to a function taking 2 arguments.

Here is my code

#include "voltra.h"

// I am trying to code your code posted above in compliance the nr3.h header file. So I have change DP to Doub,
// Vec_DP to VecDoub, and Mat_DP to MatDoub.


// I am assuming that this is g(k,t) - a user supplied function of functor that returns gk(t),
Doub g(const int k, const Doub t){
return (k == 0 ? cosh(t)+t*sin(t): 2.0*sin(t)+t*SQR(sin(t))+exp(t));
}

// I am assuming that this is ak(k,l,t,s) - a user - supplied function or functor that returns
// the (k,l) element of the matrix K(t,s)
Doub ak(const int k, const int l, const Doub t, const Doub s){
return ((k == 0) ? (l == 0 ? -exp(t+s) : -cos(t-s)) : (l == 0 ? -exp(t+s) : -t*cos(s)));
}

// Begin the main function
int main()
{
Doub pause; // dummy variable
const int N=50, M=2;
const Doub H=0.01; // Step size of 0.01;
int nn;
Doub t0=0.0; // Starting point of the integration
VecDoub g(N); // Returned corresponding abcissas
VecDoub b(N);
VecDoub t(N);

// In the voltra.h notes from NR3 the value of m is determined from the row-dimension of
// the solution matrix f and here it is specified as 2. Is this correct?

MatDoub f(M,N); // Returned solution matrix with m specified as 2 and N = 50. This will give you 2 rows and
// 50 columns. Is this correct?

voltra(t0,H,g,ak,t,f); // call to the voltra function

cout << "abcissa voltra real "
<< "voltra real" << endl;
cout << "answer1 answer2 "
<< "answer 2 answer2" << endl << endl;
cout << fixed << setprecision(6);
for(nn =0; nn<N;nn++) {
cout << setw(12) << t[nn] << setw(13) << f[0][nn];
cout << setw(13) << exp(-t[nn]) << setw(13) << f[1][nn];
cout << setw(13) << 2.0*sin(t[nn]) << endl;
}

cin >> pause; // dummy line to keep the output window visible

return 0;
}

Maybe we can work this problem out together.

Thanks
Matt

Comparing your code with mine ( i.e. NR source code ) I can spot server differences.

In order to run the code, you need several other files. These e.g. include:

nr.h
nrtypes.h
nrtypes_nr.h
nrutil.h
nrutil_nr.h

I guess you need these files to make it compile and run.

Best, "perr"

...
Below is my version of your code....

Matt: The reason it won't compile is that you declare a VecDoub named g in your main() function. Then you expect to use the function g, previously defined, as the third argument in the invocation of voltra.

So:

1. Delete or comment out the VecDoub g(N); statement.

Now it should compile, but the calculated values don't agree with the exact values because you have typos in both of your functions.

2. Change the g function to this:


Doub g(const int k, const Doub t)
{

return (k == 0 ?
cosh(t) + t * sin(t) :
2.0 * sin(t) + t * (SQR(sin(t)) + exp(t)));
}


3. Change the function ak to

Doub ak(const int k, const int l, const Doub t, const Doub s)
{
return ((k == 0) ?
(l == 0 ? -exp(t - s) : -cos(t - s)) :
(l == 0 ? -exp(t + s) : -t * cos(s)));

}


Now, at least, you should get the same answers as perr got.

In particular, the computed output for the values of parameters N and H in your code should agree closely with the exact values.

(He is apparently using NR2 C++ code and you are using NR3. I perceive that the voltra function code for the two distributions is computationally the same.)

Anyhow, use N and H values that are indicated on peer's output plots and you should see the same problems as peer stated. Namely the approximate solution diverges from the expected values of the exact solution somewhere around t=2 or t=3.

So, you can try to investigate (or think about):

Is it due to truncation error that occurs as as the trapezoidal rule is applied to step the solution or what?

For a given step size, where does the solution seem to be getting obviously worse?

What is the effect of reducing step size (and, therefore, increasing the number of steps)?

Other than ever-increasing compute time, why can the step size not be reduced indefinitely? (Roundoff error as the number of steps increases or what?)

What if a higher-order integration routine were used instead of the composite trapezoidal rule? Would an adaptive step size approach help? After all, the kernel and the solutions are both pretty well behaved, so local truncation error shouldn't be a problem, right? Or does the fact that each step of this approximate solution the integral equation actually uses all of the previous values, and this has a devastating effect on error as the previous errors accumulate?


Are there other methods that might be more suitable for solving this type of integral equation for values far from the starting point?

Maybe other, even more important, considerations could include the following, closely related questions:

Under what conditions should we expect good agreement with the known solution as the right end point value gets larger?

How could we estimate the error if we didn't know an exact solution?

What are we going to do with the results of our calculations?




Regards,

Dave

Dave, thanks. I have the routine compiling correctly now.

I have attached a graph showing my results. Yes, my results diverge.

The time step, H, is 0.01.

Other information:

Compilier - Visual C++ 2008 Express Edition
OS - Windows XP

In this example the values of M and N are defined as

const int N = 500
const int M = 2

I retrieved these values after execution of the voltra function and they came back as zero.

Why do the values of N and M change when they are typdef'd as const? I was under the impression that using the keyword const means that the values "should not" changed.

EDIT: My fault, on the values of N and M being zero. They are infact returned correctly as 500 and 2.
END EDIT:

Thanks
Matt

Hi Matt and Dave,

Now as the return value of "N" and "M" is solved, we are stuck with one final problem: the numerical solutions diverge. (C.f. my attachment in the first post and Dave's attachment "VTEX.pdf").

I contacted NR, and got reply:

"Yes, I have looked at it, but have not solved the problem about xvoltra. Last night I check the equations to make sure they were correct, and it seems that they are, so this may be an instability in the method of solution. It will take some time to investigate this."

I guess we can do nothing but wait?

Best, Per ("perr")

The equations being solved are shown in the attachment. The exact solution is
f_1(t) = exp (-t)
f_2(t) = 2 sin(t)

If you look at the second equation, you see that a term -t exp(t) from the first integral has to cancel the last term t exp(t). No matter how small you choose the step h, numerical error will accumulate to prevent an exact cancellation. So if you integrate to large enough t, the solution will become dominated by this growing exponential error. The method works fine, however, if you don't integrate too far (i.e., make N too big for a given h).

No numerical method that starts at 0 and tries to integrate to large t will work on this problem because of this instability. You would need to use a two-point boundary value type of method that incorporates the non-blowing up boundary condition at large t.

Saul Teukolsky

Hi Saul,

Thank you for your reply.

As I understand your post, there is nothing wrong with the NR method "voltra.cpp". The divergence is due to a cancellation that numerically only almost cancel, and hence the solution will become dominated by this growing exponential error. (I then believe this spesific example in "xvoltra.cpp" is a poor choice).

So I can safely apply "voltra.cpp" for my realife physical system?

Best regards, Per ("perr")

That's right.

Saul