Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Computing Integration of Bessel Function

  1. Mar 26, 2012 #1
    I tried to compute this exact solution, but faced difficulty if the value of [itex] η[/itex] approaching to [itex] ζ [/itex]. Let say the value of [itex] ζ [/itex] is fix at 0.5 and the collocation points for η is from 0 to 1.

    [tex] θ(η,ζ)=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(ζ-η)[/tex]

    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
     
  2. jcsd
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted



Similar Discussions: Computing Integration of Bessel Function
Loading...