clenshawm.m

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

%M阶（M>1,M是整数)MarcumQ函数计算


clc;
close all;
clear all;

M=85;
Belta=10;
Alpha=[1:1:20];
AAlpha=Alpha;
L=length(Alpha);
marcumq1=zeros(1,L);   %直接用marcumq函数计算，作为参考
q1=zeros(1,L);         %用三项求和法
epsilon=zeros(1,L);   %相对误差
n=0;
for Alpha=1:1:20
n=n+1;
%先处理Alpha>Belta的情况
if(Alpha>Belta)
%用alpha表示两者中较小者，用belta表示两者中较大者   
alpha=Belta;
belta=Alpha;
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:M
    A=-2*k/(belta*belta);
    B=(alpha/belta)^2;
    y1(k)=A*y1(k+1)+B*y1(k+2)+1;
end

%按k递增的顺序计算
if y1(M+1)*y1(M+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);
FM_after=caita^(M+1)*besseli(M+1,alpha*belta);
FM=caita^M*besseli(M,alpha*belta);
%级数求和的结果
f_M2=FN-B*FN_1*y1_2(N1-1+3)-FN*y1_2(N1-2+3)+FM_after*y1_2(M-1+3)+(y1_2(M-2+3)-(-2*M/(belta*belta))*y1_2(M-1+3))*FM;
    
end



if y1(M+1)*y1(M+2)>=0
FM_after=caita^(M+1)*besseli(M+1,alpha*belta);
FM=caita^M*besseli(M,alpha*belta);
f_M1=B*y1(M+2)*FM+y1(M+1)*FM_after+1*FM;
end  

%结果及误差
if(y1(M+1)*y1(M+2)>=0)
q1(n)=1-exp(-0.5*(alpha*alpha+belta*belta))*f_M1;
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(M+1)*y1(M+2)<0)
q1(n)=1-exp(-0.5*(alpha*alpha+belta*belta))*f_M2;
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)


if (Belta>Alpha)
    %用（4.37）或（4.38）计算
    %用alpha表示两者中较小者，用belta表示两者中较大者
alpha=Alpha;
belta=Belta;
caita=alpha/belta;


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

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

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


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


if(y1(1)*y1(2)>=0)
for i=1-M:-1
    f_1=f_1+caita^i*besseli(-i,alpha*belta);
end
q1(n)=exp(-0.5*(alpha*alpha+belta*belta))*f_1;
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)
for i=1-M:-1
    f_2=f_2+caita^i*besseli(-i,alpha*belta);
end
q1(n)=exp(-0.5*(alpha*alpha+belta*belta))*f_2;
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';