0

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;
}
3
Contributors
3
Replies
4
Views
6 Years
Discussion Span
Last Post by sniper6560
0

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

Have something to contribute to this discussion? Please be thoughtful, detailed and courteous, and be sure to adhere to our posting rules.