std::cyl_bessel_i in C++17

In Mathematics, the differential equation,

x^{2}\frac{d^{2}y}{dx^{2}} + x\frac{dy}{dx} + (x^{2} - v^{2})y = 0

is of particular importance. The solution of the equation is a function of the parameter v. For integral and half-integral values of v, the solutions of particular interest and are called as Cylindrical Bessel’s Functions, named after the famous German mathematician, Friedrich Wilhelm Bessel. The reason for the requirement of v to be integral or half-integral will be clear in the explanation given below.



Since this is a second-order differential equation, there must be two linearly independent solutions, called as the first kind and the second kind. Thus, the differential equation can be solved by hand using the Frobenius Method. Bessel functions of the first kind, for complex arguments x, are called Modified Bessel Functions of the first kind and are denoted by I_{v}(x). The application of the method produces an infinite series containing terms of x and v, given by,

 I_{v}(x) = \sum_{k = 0}^{\infty} \frac{1}{k!\Gamma (k+v+1)}\left(\frac{x}{2}\right)^{2k+v}

Since the expression contains Gamma Function, which can only be calculated for integral and half-integral values, thus the parameter v must be integral or half-integral.

C++17 (GCC 7.1) Standard Library cmath gives functions which compute the value of the Cylindrical Bessel’s Function of the first kind (std::cyl_bessel_j) (not discussed here, but is very similar to what we have discussed) and the value of Regular Modified Bessel Functions (std::cyl_bessel_i). Both have appreciable accuracy for small inputs and can be used in various engineering applications.

Examples:

Input: x = 2.798465, v = 0
Output: 4.152234090041574

Input: x = 3.04513, v = 0.5
Output: 4.792979684692604

Note: The following source code should only be run on C++17 and above. The running sample of the given code can be checked here. To run for a different input, please visit the link and click “Edit” in the bottom right corner.

filter_none

edit
close

play_arrow

link
brightness_4
code

// C++17 code for bessel function
#include <bits/stdc++.h>
using namespace std;
  
// Compute the answer from the formulae for first 10 terms
long double answer(long double x, long double v)
{
  
    long double ans_by_expansion = 0;
    long double fact = 1;
  
    for (int k = 0; k < 10; fact = fact * (++k)) {
        ans_by_expansion += pow((x / 2), (2 * k)) / pow(fact, 2);
        cout << "ans_by_expansion till term k = ";
        cout << k << " is " << ans_by_expansion << "\n";
    }
  
  return ans_by_expansion;
}
  
// Driver code
int main()
{
    long double x = 2.798465;
    long double v = 0;
  
    // Compute the Regular Modified Bessel Function
    // for v = 0, x = 2.798465
    long double ans_by_function = cyl_bessel_i(v, x);
  
    cout << setprecision(15) << fixed;
    cout << "The answer by function for "
         << "Regular_Modified_Bessel_Function" << endl
         << "(" << v << ", " << x << ") = "
         << ans_by_function << "\n";
  
    // calculate answer by expansion
    long double ans_by_expansion = answer(x, v);
  
    cout << "Absolute Error in answer by both the methods is = ";
    cout << abs(ans_by_expansion - ans_by_function) << "\n";
  
    return 0;
}

chevron_right


Output:

The answer by function for Regular_Modified_Bessel_Function
(0.000000000000000, 2.798465000000000) = 4.152234090041574
ans_by_expansion till term k = 0 is 1.000000000000000
ans_by_expansion till term k = 1 is 2.957851589056250
ans_by_expansion till term k = 2 is 3.916147300248771
ans_by_expansion till term k = 3 is 4.124614053687001
ans_by_expansion till term k = 4 is 4.150123238967278
ans_by_expansion till term k = 5 is 4.152120966924739
ans_by_expansion till term k = 6 is 4.152229612892962
ans_by_expansion till term k = 7 is 4.152233953968095
ans_by_expansion till term k = 8 is 4.152234086767796
ans_by_expansion till term k = 9 is 4.152234089977698
Absolute Error in answer by both the methods is = 0.000000000063876


My Personal Notes arrow_drop_up

Check out this Author's contributed articles.

If you like GeeksforGeeks and would like to contribute, you can also write an article using contribute.geeksforgeeks.org or mail your article to contribute@geeksforgeeks.org. See your article appearing on the GeeksforGeeks main page and help other Geeks.

Please Improve this article if you find anything incorrect by clicking on the "Improve Article" button below.


Article Tags :
Practice Tags :


Be the First to upvote.


Please write to us at contribute@geeksforgeeks.org to report any issue with the above content.