Computing Integration of Bessel Function

[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

 

[jvicens]

Thank you for sharing your experience and code with us. It seems like you have put a lot of effort into computing this exact solution and I appreciate your dedication.

From my understanding, the main issue you are facing is that you are not getting the expected value of θ when η = 0.5. This could be due to the singularity at that point, as you have mentioned. One possible solution could be to use a different integration method that is more suitable for handling singularities, such as the adaptive quadrature method.

Additionally, it might be helpful to check your code for any errors or bugs, as even a small mistake can greatly affect your results. You can also try to break down your code into smaller parts and check each step to see if it is producing the expected results.

I hope these suggestions will help you in computing the correct value of θ at η = 0.5. Keep up the good work and don't hesitate to reach out if you need further assistance.