# Computing Integration of Bessel Function

1. Mar 26, 2012

### sugaku

I tried to compute this exact solution, but faced difficulty if the value of $η$ approaching to $ζ$. Let say the value of $ζ$ is fix at 0.5 and the collocation points for η is from 0 to 1.

$$θ(η,ζ)=e^{-ε\frac{η}{2}} \left\{ e^{-η}+\left(1-\frac{ε^2}{4}\right)^{1/2} η \int_η^ζ e^{-τ}\frac{I_1 \left\{\left[ \left(1-\frac{ε^2}{4}\right)\left(τ^2-η^2\right) \right]^{1/2} \right\}} {\left(τ^2-η^2\right)} \right\} U(ζ-η)$$

These are the values that is suppose to appear, but only when η=0.5 θ=0.295778, i don't manage to get that value, others is ok. I used trapz command in matlab to calculate the area.

η=0.0 θ=1.000000
η=0.1 θ=0.915287
η=0.2 θ=0.831763
η=0.3 θ=0.749758
η=0.4 θ=0.669587
η=0.5 θ=0.295778
η=0.6 θ=0.000000
η=0.7 θ=0.000000
η=0.8 θ=0.000000
η=0.9 θ=0.000000
η=1.0 θ=0.000000

I do suspect that the integration of Bessel function is not simply become 0 when η=0.5 (approach to singularity to that point). I do appreciate if someone could give some advice.

Here I attach the matlab program that I wrote. Thank you in advance

format short
%analytic solution
tic

ita=0:0.1:1; m=11;

ep=0.1;
zeta=0.5;

area=zeros(1,m);
%kira integration dahulu
for i=1:m
if ita(i)<=zeta
tau=linspace(ita(i)+0.000001,0.5,100000);
%argument for bessel function
a=(1-(ep^2)/4);b=(tau.^2-ita(i)^2);
Z=(a*b).^(1/2);
%Modified bessel function
func=@(tau) (exp(-tau).*besseli(1,Z))./sqrt(b);
area(i)=trapz(tau,func(tau));
else
area(i)=0;
end
end

Theta=zeros(1,m);
for i=1:m
if ita(i)<=zeta
Theta(i)=exp((-ita(i)./2)*ep)*(exp(-ita(i))+sqrt(a)*ita(i).*area(i));
else
Theta(i)=0;
end
end

plot(ita,Theta);
axis([0 2.2 0 1]);

tableresult(:,1)=ita';
tableresult(:,2)=Theta';

disp(' x Analytic')
disp('')
disp(tableresult);
toc