Trying of implementation of Bessel function of first kind

This is my first article about mathematics in dlang written in English. Maybe, I made several mistakes in this short article because I only began to learn English.

So, I will to narrate about how I wrote Bessel function of first kind in D programming language.

If you want to read about this — welcome under shortcut ;)

Once, I tried to program Bessel function of first kind, in that time, when I used Icon language for solving any my mathematical problems. Then I did this function without any understanding how she works. Unfortunately, I couldn’t understand this function because I was too weak in special functions…

Let’s sort out with Bessel function of first kind.

Bessel functions are the class of function that are canonical solutions of differential equation, which is defined as follows:

where variable alpha is arbitrary real number, called order of Bessel function. Mostly use integer value of order of these function.

For Bessel functions, exists Taylor series that I used for own implementation of these functions. It’s looks like this:

where variable x is argument for Bessel function, variable alpha — order of Bessel first kind function and Г is Euler’s gamma function.

D has gamma function in his own standart library (She can be found in module std.mathspecial). We will to use her and for getting arbitrary precision of Bessel function we will apply variables with D’s type real.

However, for getting complete implementation of function need to implement factorial function, which can be implemented by iteration like this:

real factorial(real N)
{
	auto P = 1;

	for (real i = 1; i < N; i++)
	{
		P *= i;
	}

	return P;
}

Now, we can go to Bessel function, which looks like this:

real bessel_1(real x, real alpha)
{
	const real multiplicationFactor = (0.5 * x) ^^ alpha;
	real S = 0.0;

	auto numerator(real m)
	{
		return ((-1.0) ^^ m) * ((0.5 * x) ^^ (2 * m));
	}

	auto denominator(real m)
	{
		return factorial(m) * gamma(m + alpha + 1);
	}

	for (real m = 0.0; m < 20; m++)
	{
		S += numerator(m) / denominator(m);
	}

	return S * multiplicationFactor;
}

I simplified Taylor series by extracting factor (x / 2) ^^ alpha, because he can be defined as constant for concrete alpha. Also, I partition full expression for Bessel function on numerator and denominator, which are defined separately as local function inside bessel_1 procedure.

For sufficient precision is not necessary to compute infinite Taylor series. Therefore, we can adjust precision by selecting number of repetition of cycle for (For me is sufficient that this number equals to 20). You can self to select precision empirically.

After computing sum need to multiply her by our constant, which we called multiplicationFactor. This provides full value for Bessel function of first kind for any arbitrary order.

Now we can to test her:

void main()
{
	writefln("%.20f", bessel_1(1.5, 4.0));
}

In this short example, I used format specifier «%.20» for getting significand with 20 digits after point, of course, you can use any format specifier for your mathematical task.

Full code of example:

import std.mathspecial;
import std.stdio;


real factorial(real N)
{
	auto P = 1;

	for (real i = 1; i < N; i++)
	{
		P *= i;
	}

	return P;
}


real bessel_1(real x, real alpha)
{
	const real multiplicationFactor = (0.5 * x) ^^ alpha;
	real S = 0.0;

	auto numerator(real m)
	{
		return ((-1.0) ^^ m) * ((0.5 * x) ^^ (2 * m));
	}

	auto denominator(real m)
	{
		return factorial(m) * gamma(m + alpha + 1);
	}

	for (real m = 0.0; m < 20; m++)
	{
		S += numerator(m) / denominator(m);
	}

	return S * multiplicationFactor;
}

void main()
{
	writefln("%.20f", bessel_1(1.5, 4.0));
}

It’s all for today, I hope that that article will to be not last.

P.S: The author thanks to Balakina Daria for helping with the correction of grammatical mistakes in the article.