Anterior

Funciones de Bessel

Raíces de una ecuación

Procedimiento del punto medio

Máximos y mínimos de difracción de auna abertura circular

public abstract class Raiz {
	private final double CERO=1e-10;
	private final double ERROR=0.0001;
	private final int MAXITER=100;
	protected final int MAXRAICES=20;
	protected double raices[]=new double[MAXRAICES];
	protected int iRaiz=0;

protected double puntoMedio(double a, double b)throws RaizExcepcion{
	double m, ym;
	int iter=0;
	do{
		m=(a+b)/2;
		ym=f(m);
		if(Math.abs(ym)<CERO) break;
		if(Math.abs((a-b)/m)<ERROR) break;
		if((f(a)*ym)<0) b=m;
		else a=m;
		iter++;
	}while(iter<MAXITER);
	if(iter==MAXITER){
		throw new RaizExcepcion("No se ha encontrado una raíz");
	}
	return m;
}

private void explorar(double xIni, double xFin, double dx){
	double y1, y2;
	iRaiz=0;
	y1=f(xIni);
	for(double x=xIni; x<xFin; x+=dx){
		y2=f(x+dx);
//Uno de los extremos del intervalo es raíz
		if(Math.abs(y1)<CERO && iRaiz<MAXRAICES){
			raices[iRaiz++]=x;
			y1=y2;
			continue;
		}
//no hay raíz en este intervalo
		if(y1*y2>=0.0){
			y1=y2;
			continue;
		}
//hay una raíz en este intervalo
		if(iRaiz<MAXRAICES){
			try{
				raices[iRaiz]=puntoMedio(x, x+dx);
				iRaiz++;
			}catch(RaizExcepcion ex){
				System.out.println(ex.getMessage());
			}
		}
		y1=y2;
	}
}

public double[] hallarRaices(double ini, double fin, double paso){
	explorar(ini, fin, paso);
	double solucion[]=new double[iRaiz];
	for(int i=0; i<iRaiz; i++){
		solucion[i]=(double)Math.round(raices[i]*1000)/1000;
	}
	return solucion;
}

abstract protected double f(double x);
}

class RaizExcepcion extends Exception {
	public RaizExcepcion(String s) {
	super(s);
}
}
public class Maximos extends Raiz {
protected double f(double x){
	return Bessel.bessj(2, x);
}

}
public class Minimos extends Raiz {
protected double f(double x){
	return Bessel.bessj1(x);
}


public class Difraccion {

public static void main(String[] args) {
	double[] raices=new Maximos().hallarRaices(0.0, 30.0, 1.0);
	System.out.print("Máximos");
	for(int i=0; i<raices.length; i++){
		System.out.print(" "+raices[i]);
	}
	System.out.println("");
	raices=new Minimos().hallarRaices(0.0, 30.0, 1.0);
	System.out.print("Mínimos");
	for(int i=0; i<raices.length; i++){
		System.out.print(" "+raices[i]);
	}
}
}
public class Bessel {
public static double bessj0(double x){
//Returns the Bessel function J0(x) for any real x.
	double ax,z;
	double xx,y,ans,ans1,ans2; //Accumulate polynomials in double precision.
	if ((ax=Math.abs(x)) < 8.0) { //Direct rational function t.
		y=x*x;
		ans1=57568490574.0+y*(-13362590354.0+y*(651619640.7+y*(
-11214424.18+y*(77392.33017+y*(-184.9052456))))); ans2=57568490411.0+y*(1029532985.0+
y*(9494680.718+y*(59272.64853+y*(267.8532712+y*1.0)))); ans=ans1/ans2; } else { //Fitting function (6.5.9). z=8.0/ax; y=z*z; xx=ax-0.785398164; ans1=1.0+y*(-0.1098628627e-2+y*(0.2734510407e-4 +
y*(-0.2073370639e-5+y*0.2093887211e-6))); ans2 = -0.1562499995e-1+y*(0.1430488765e-3 +y*(-0.6911147651e-5+ y*(0.7621095161e-6-y*0.934935152e-7))); ans=Math.sqrt(0.636619772/ax)*(Math.cos(xx)*ans1-z*Math.sin(xx)*ans2); } return ans; } } public static double bessj1(double x){ //Returns the Bessel function J1(x) for any real x. double ax,z; double xx,y,ans,ans1,ans2; //Accumulate polynomials in double precision. if ((ax=Math.abs(x)) < 8.0) { //Direct rational approximation. y=x*x; ans1=x*(72362614232.0+y*(-7895059235.0+y*(242396853.1+y*(-2972611.439+ y*(15704.48260+y*(-30.16036606)))))); ans2=144725228442.0+y*(2300535178.0
+y*(18583304.74+y*(99447.43394+y*(376.9991397+y*1.0)))); ans=ans1/ans2; } else { //Fitting function (6.5.9). z=8.0/ax; y=z*z; xx=ax-2.356194491; ans1=1.0+y*(0.183105e-2+y*(-0.3516396496e-4
+y*(0.2457520174e-5+y*(-0.240337019e-6)))); ans2=0.04687499995+y*(-0.2002690873e-3+y*(0.8449199096e-5+ y*(-0.88228987e-6+y*0.105787412e-6))); ans=Math.sqrt(0.636619772/ax)*(Math.cos(xx)*ans1-z*Math.sin(xx)*ans2); if (x < 0.0) ans = -ans; } return ans; } public static double bessj(int n, double x){ //Returns the Bessel function Jn(x) for any real x and n > 2. final double ACC=40.0; // Make larger to increase accuracy. final double BIGNO=1.0e10; final double BIGNI=1.0e-10; //int j,jsum,m; boolean jsum; int j, m; double ax,bj,bjm,bjp,sum,tox,ans; if(n==0) return bessj0(x); if(n==1) return bessj1(x); if(n==-1) return -bessj1(x); // if (n < 2) nrerror("Index n less than 2 in bessj"); ax=Math.abs(x); if (ax == 0.0) return 0.0; else if (ax > (double) n) { //Upwards recurrence from J0 and J1. tox=2.0/ax; bjm=bessj0(ax); bj=bessj1(ax); for (j=1;j<n;j++) { bjp=j*tox*bj-bjm; bjm=bj; bj=bjp; } ans=bj; } else { //Downwards recurrence from an even m here computed. tox=2.0/ax; m=2*((n+(int) Math.sqrt(ACC*n))/2); jsum=false; //jsum will alternate between 0 and 1;
//when it is 1, we accumulate in sum the even terms in (5.5.16). bjp=ans=sum=0.0; bj=1.0; for (j=m;j>0;j--) {// The downward recurrence. bjm=j*tox*bj-bjp; bjp=bj; bj=bjm; if (Math.abs(bj) > BIGNO) { //Renormalize to prevent overflows. bj *= BIGNI; bjp *= BIGNI; ans *= BIGNI; sum *= BIGNI; } if (jsum) sum += bj; // Accumulate the sum. jsum=!jsum; //Change 0 to 1 or vice versa. if (j == n) ans=bjp; //Save the unnormalized answer. } sum=2.0*sum-bj; //Compute (5.5.16) ans /= sum; //and use it to normalize the answer. } boolean bRes=((n&1)==0)? false:true; return (x < 0.0) && bRes ? -ans : ans; } }

Referencias

Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P. Numerical Recipes in C, Second edition,  Special functions. Bessel functions of integer order  Chapter 6º. pp. 230. Cambridge University Press. Código en C adaptado por el autor al lenguaje Java

Anterior