bessel function 2nd kind Yn(x)

sniper6560 0 Tallied Votes 1K Views Share

this code is for bessel function Y n (x) 2nd kind I'm in the level 4 at computer
engineering my prof. gave me that assignment , I want help to improve it
to get an accurate results thanks!

#include<iostream>//bessel function Y n (x) 2nd kind  
using namespace std;
#include<cmath>
long double Segma1 (long double,long double);
long double Segmafact (long double);
long double Segma2 (long double,long double);
long double HMMM (long double);
long double factorials (long double,long double);
long double Jn (long double,long double);
int main()
{
	long double m=2,c,s,p,b,bessel1,segma1,segma2,Yn;
	const long double pi=3.14159265358,EulerConstant=0.577215664901532860606512;
		for(double n=0;n<=5;n++)
{
	cout<<"n= "<<n<<endl<<endl;
	for(double x=0;x<=20;x+=0.5)
	{
	bessel1= Jn(x,n);
	segma1=Segma1(n,x);
	segma2=Segma2(n,x);
	c=log(x/2)+EulerConstant;
	s=m/pi;
	p=(pow(x,n))/pi;
	b=(pow(x,(-n)))/pi;
	Yn=(s*bessel1*c)+(p*segma1)-(b*segma2);
	cout<<"Yn("<<x<<") =  "<<Yn<<endl;
	Yn=0;
	}
		}
	return 0;
}
long double Jn (long double x,long double n)
{
	long double k=0,t,z,b,r,m=0,E=0.000001;
	
	for( ; ; ){
		t=n+k;
        r=n+(2*k);
        z=factorials(k,t);
		b=z*(pow(-1,k))*(pow((x/2),r));
		
		if(fabs(b)<E)
			break;
		else
			m+=b;
			k++;
	}
	return m;
}
long double Segma1 (long double n,long double x)
{
	long double f,g,Hm,Hh,h,Segma1A=0,o=-1,p=2,m=0;
	for( ;  ;)
	{
		h=(m+n);
		f=Segmafact(h);
	    g=Segmafact(m);
		Hm=HMMM(m);
		Hh=HMMM(h);
		Segma1A+=(pow(o,(m-1))*(Hm+Hh)*pow(x,(2*m)))/(pow(p,(2*m+n))*g*f);
		if(m==99)
			break;
		m++;
	}
return Segma1A;
}
long double Segma2 (long double n,long double x)
{
	long double t,f,g ,Segma2B=0;
	for(long double m=0;m<=(n-1);m++)
	{
		t=n-m-1;
		f=Segmafact(t);
	    g=Segmafact(m);
		Segma2B+=(f*pow(x,(2*m)))/(pow(2,(2*m-n))*g);
	}
	return Segma2B;
}
long double factorials (long double k,long double t)
{
	long double fact1=1, fact2=1, last;
	for(long double j=1;j<=k;++j){
if(k==0)
fact1=1;
else
{fact1*=j;}}
	 for (long double i=1;i<=t; ++i)
	 {
	  fact2*=i;
	 }
	 last=1/(fact1*fact2);
	 return last;
}
long double Segmafact (long double n){
	long double fact=1;
	if(n!=0)
	{
		if(n>0){
	  for(long double j=1; j<=n; j++)
			{
				fact*=j;
			}
		}
		else
			{
				for(long double j=n; j<=(fabs(n)); j++)
			{
				fact*=j;
			}
				fact*=-1;
			}
	}
	else
		fact=1;	
		return fact;
	} 
long double HMMM (long double s)
{   
	long double o=1, l,h=0;
	for(long double n=1;n<=s;n++){
		l=(o/n);
		h+=l;
	}
	return h;
}
daviddoria 334 Posting Virtuoso Featured Poster
pradev 0 Newbie Poster

Hello,

As far as I know, at the origin Yn(x) is singular.
Consequently, the result you get for x=0 is wrong anyway.

There is a recursion relation for J, n/xJn=1/2(Jn-1 + Jn+1), so once
you have calculated two Jn=0,1 you can reuse these to calculate the Jn n>1

Use M_PI instead of your own definition.

Hope this helps.

Johan

sniper6560 0 Newbie Poster

thank you for that Johan !

Be a part of the DaniWeb community

We're a friendly, industry-focused community of developers, IT pros, digital marketers, and technology enthusiasts meeting, networking, learning, and sharing knowledge.