Difusión. Ley de Fick

Procedimiento numérico

Funciones de Bessel. Integral de Simpson

$\frac{n (x, τ)}{n_{0}} = \frac{1}{2 τ} \exp (\frac{- x^{2}}{4 τ}) \int_{0}^{1} \exp (\frac{- ζ^{2}}{4 τ}) \cdot I_{0} (\frac{x \cdot ζ}{2 τ}) \cdot ζ \cdot d ζ$

//método de Simpson

public double integral(double a, double b, int n){
	if(n%2==1) n++;
	double h=(b-a)/n;
	double suma=f(a)+f(b);
	for(int i=1; i<n; i+=2){
		suma+=4*f(a+i*h);
	}
	for(int i=2; i<n; i+=2){
		suma+=2*f(a+i*h);
	}
	return (suma*h/3);
}

double f(double z){
	double temp=Math.exp(-z*z/(4*t))*bessi0(x*z/(2*t))*z;
	return temp;
}


public double bessi0(double x) {
//Returns the modifed Bessel function I0(x) for any real x.
	double ax,ans;
	double y; //Accumulate polynomials in double precision.
	if ((ax=Math.abs(x)) < 3.75) { //Polynomial t.
		y=x/3.75;
		y*=y;
		ans=1.0+y*(3.5156229+y*(3.0899424+y*(1.2067492
+y*(0.2659732+y*(0.360768e-1+y*0.45813e-2)))));
	} else {
		y=3.75/ax;
		ans=(Math.exp(ax)/Math.sqrt(ax))*(0.3989422
8+y*(0.1328592e-1+y*(0.225319e-2+
		y*(-0.157565e-2+y*(0.916281e-2+y*(-0.2057706e-1+y*(0.2635537e-1
		+y*(-0.1647633e-1+y*0.392377e-2))))))));
	}
	return ans;
}

//se representa la función

 y=Math.exp(-x*x/(4*t))*integral(0.0, 1.0, 40)/(2*t);

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 Jav