- #1
sugaku
- 17
- 0
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
[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