clenshaw.asv

用“三项求和法”计算正整数阶MarcumQ函数

%一阶MarcumQ函数计算
clc;
close all;
clear all;

M=1;
Belta=17.6022;
Alpha=[0.5:0.5:20];
AAlpha=Alpha;
L=length(Alpha);
marcumq1=zeros(1,L);   %直接用marcumq函数计算，作为参考
q1=zeros(1,L);         %用三项求和法
epsilon=zeros(1,L);    %相对误差
n=0;
for Alpha=0.5:0.5:20
n=n+1;
%用alpha表示两者中较小者，用belta表示两者中较大者
alpha=Alpha;
belta=Belta;
if(Alpha>Belta)
alpha=Belta;
belta=Alpha;
end

caita=alpha/belta;

%确定N
flag=1;
i=100;
while (i<inf&&flag==1)
    term=caita^i*besseli(i,alpha*belta);
    if(term<10^-100)
        N1=i;
        flag=0;
        break;
    end
    i=i+100;
end

%计算yk
y1=zeros(1,N1+2);
y1(N1+2)=0;
y1(N1+1)=0;
for k=N1:-1:1
    A=-2*k/(belta*belta);
    B=(alpha/belta)^2;
    y1(k)=A*y1(k+1)+B*y1(k+2)+1;
end

%按k递增的方向重新计算yk
if y1(1)*y1(2)<0
    y1_2=zeros(1,N1+2);
y1_2(1)=0; %y-2
y1_2(2)=0;  %y-1
B=(alpha/belta)^2;
for k=3:1:N1+2
    A=-2*(k-3)/(belta*belta);
    y1_2(k)=1/B*(y1_2(k-2)-A*y1_2(k-1)-1);
end

FN=(alpha/belta)^N1*besseli(N1,alpha*belta);
FN_1=(alpha/belta)^(N1-1)*besseli(N1-1,alpha*belta);
f_2=FN-B*FN_1*y1_2(N1-1+3)-FN*y1_2(N1-2+3);   %三项求级数和
    
end


if( y1(1)*y1(2)>=0)   
F0=besseli(0,alpha*belta);
F1=(alpha/belta)*besseli(1,alpha*belta);
f_1=B*F0*y1(2)+F1*y1(1)+F0*1;
end



%结果及误差，与Alpha和Belta的相对大小有关
if (Belta>Alpha)
    %用（4.37）或（4.38）计算
if(y1(1)*y1(2)>=0)

%fprintf('方法一相对误差为%e\n',error1_1);
q1(n)=exp(-0.5*(alpha*alpha+belta*belta))*f_1;
%disp('方法一计算MarcumQ函数结果为');
%q1(Alpha)
%disp('用marcumq函数计算MarcumQ函数结果为');
marcumq1(n)=marcumq(Alpha,Belta,M);
epsilon(n)=(marcumq1(n)-q1(n))/marcumq1(n);
fprintf('当Alpha=%e时，方法一最终结果的相对误差为%e\n',Alpha,epsilon(n));
end

if(y1(1)*y1(2)<0)

%fprintf('方法二相对误差为%e\n',error1_2);
q1(n)=exp(-0.5*(alpha*alpha+belta*belta))*f_2;
%disp('方法二计算MarcumQ函数结果为');
%q1_2
%disp('用marcumq函数计算MarcumQ函数结果为');
marcumq1(n)=marcumq(Alpha,Belta,M);
epsilon(n)=(marcumq1(n)-q1(n))/marcumq1(n);
fprintf('当Alpha=%e时，方法二最终结果的相对误差为%e\n',Alpha,epsilon(n));
end

end

if(Belta<Alpha)
    %用式（4.36）或(4.39)计算
    if(y1(1)*y1(2)>=0)

%fprintf('方法一相对误差为%e\n',error1_1);
q1(n)=1-exp(-0.5*(alpha*alpha+belta*belta))*(f_1-besseli(0,alpha*belta));
%disp('方法一计算MarcumQ函数结果为');
%q1(Alpha)
%disp('用marcumq函数计算MarcumQ函数结果为');
marcumq1(n)=marcumq(Alpha,Belta,M);
epsilon(n)=(marcumq1(n)-q1(n))/marcumq1(n);
fprintf('当Alpha=%e时，方法一最终结果的相对误差为%e\n',Alpha,epsilon(n));
end

if(y1(1)*y1(2)<0)

%fprintf('方法二相对误差为%e\n',error1_2);
q1(n)=1-exp(-0.5*(alpha*alpha+belta*belta))*(f_2-besseli(0,alpha*belta));
%disp('方法二计算MarcumQ函数结果为');
%q1(Alpha)
%disp('用marcumq函数计算MarcumQ函数结果为');
marcumq1(n)=marcumq(Alpha,Belta,M);
epsilon(n)=(marcumq1(n)-q1(n))/marcumq1(n);
fprintf('当Alpha=%e时，方法二最终结果的相对误差为%e\n',Alpha,epsilon(n));
end


end

end  %end for


 marcumq1=marcumq1';
 AAlpha=AAlpha';
 q1=q1';
 epsilon=epsilon';