Hi
I have a very stupid problem. I would like to use the modified Bessel function with fractional order. Therefore I use bessik() from NR. First I compared the result for order 0 and argument x with bessi0(x). I expected the same result. But the results are different. Do I use the wrong function? But this is the modified bessel function... I am puzzled
thanks for any help
best regards
florian

I would like to use the modified Bessel function with fractional order. Therefore I use bessik() from NR. First I compared the result for order 0 and argument x with bessi0(x).


#include <iostream>
#include <iomanip>
#include "nr.h"

using namespace std;

int main()
{
DP x, xnu, xri, xrk, xrip, xrkp;
DP result;

x = 3.5;
xnu = 0.0;

result = NR::bessi0(x);

cout << " bessi0("
<< fixed << setprecision(1)
<< x << ") = "
<< scientific << setprecision(6)
<< result << endl;

NR::bessik(x, xnu, xri, xrk, xrip, xrkp);

cout << "From bessik("
<< fixed << setprecision(1)
<< x << ", " << xnu << "), ri = "
<< scientific << setprecision(6)
<< xri << endl;


return 0;
}

Output on my Linux system (g++ version 4.1.2):


bessi0(3.5) = 7.378203e+00
From bessik(3.5, 0.0), ri = 7.378203e+00


Regards,

Dave

thanks for the quick reply. Ok obviously there is a typo in my code. I used the original NR stuff, but I copied just all necessary functions in one. I attached the file. Maybe I use an old function?
regards
florian

[editor's note: it is not permitted to post NR source code on the forum]

thanks for the quick reply. Ok obviously there is a typo in my code. I used the original NR stuff, but I copied just all necessary functions in one. I attached the file. Maybe I use an old function?
regards
florian

It has been said many, many (many) times (by Bill Press, Saul Teukolsky, and, I'm sure, many others):


Most reports of user problems are from typing errors. I don't have enough experience over enough time with enough users to justify the "most" part of it, but I will accept that statement with the minor amplification: "or copy/paste errors."

Somewhere around line 56 of your file:


ril=FPMIN;
ripl=h*ril; //<--- this statement was missing
ril1=ril;
rip1=ripl;



General comments:

1. You are apparently using the C functions in a C++ program. My suggestion is that if you are going to use C functions, write a C program. If you are going to write a C++ program, then use C++ functions. That way, you won't have to copy anything or edit anything. Just use the functions as supplied unless and until you really (really) need to change something for your particular application. Additionally, I add the opinion that commenting out validity tests in NR-supplied functions seems particularly counter-productive in trying to produce a generally useful program.

2. When you ask a question, I think that you increase the chances of getting a helpful response with fewest iterations if you tell us a few specific things.

[begin sermon]
Things like:

What, exactly did you do. (Stuff like, "Compiled the C functions all together as a C++ program.")

Show us a small main() function so that we can see how you tested the NR function(s). How do we know that you used the functions correctly? (Especially given that you commented out the places where the functions checked for valid input values.) We might offer suggestions based on guesses, but why should we have to guess?

Tell us exactly what the program told you. (Stuff like, "It gave the correct answer for bessi0() result but not for the bessik()". I mean, surely you wouldn't test the functions without knowing what the correct answers were, right? Anyhow, show us exactly what your version of your test program produced as output.

You might also tell us what compiler you used. Sometimes it makes a difference to people who would like to help.
[/end sermon]


After adding the missing line to your code for bessik, and changing my previous version of the main() function to the following:


#include <iostream>
#include <iomanip>
#include "nr.h"

using namespace std;

int main()
{
float x, xnu, xri, xrk, xrip, xrkp;
float result;

x = 3.5;
xnu = 0.0;

result = bessi0(x);

cout << " bessi0("
<< fixed << setprecision(1)
<< x << ") = "
<< scientific << setprecision(6)
<< result << endl;

bessik(x, xnu, &xri, &xrk, &xrip, &xrkp);

cout << "From bessik("
<< fixed << setprecision(1)
<< x << ", " << xnu << "), ri = "
<< scientific << setprecision(6)
<< xri << endl;


return 0;
}



I compiled this along with your file, and here's what I get (on my Linux system with g++ version 4.1.2):


bessi0(3.5) = 7.378203e+00
Bessi is runing
From bessik(3.5, 0.0), ri = 7.378203e+00


(You misspelled "running," but I won't deduct any points for that. This time. Also, I will offer the observation that having a statement in a function named bessik that prints out misleading information---the name of another function---could be a little confusing in the grand scheme of things. Details, details...)

Regards,

Dave