关于运算精度问题！

跑出来不对，wj1,wj2是1*1的维度！

应该是少了冒号我加上去看看！

修改了你的代码，但是跑出来的解是不对的，两个解大于零，在opw的左右，这个最好用二分法解，fsolve解两解的方程会出错，我在尝试二分法!

乱云而已 发表于 2016-5-15 16:37
跑出来不对，wj1,wj2是1*1的维度！

110.665967971709        164504.634259209
386.787331398439        65793.9771025985
514.759935758331        47564.6630943982
600.775343131164        38762.4458188561
665.664403517200        33380.7741209675
717.751912872416        29677.5069148783
761.242049961360        26938.9477739223
798.559528367342        24812.8053857315
831.230997565907        23103.0828437254
860.280040290540        21691.1718697730

前十组解, 你看看对不对

乱云而已 发表于 2016-5-15 16:52
修改了你的代码，但是跑出来的解是不对的，两个解大于零，在opw的左右，这个最好用二分法解，fsolve解两解的 ...

function [w, wj1, wj2] = Untitled
m=500;n=5;rou=1;alpha=0.0023;
s=rou*m*(n-1);
ucl=chi2inv(1-alpha/2,n-1);
lcl=chi2inv(alpha/2,n-1);
opw=rou*m*(n-1)^2*((log(ucl/lcl))/(ucl-lcl));
for t=1.1:541.1
    f = @(w)1-chi2cdf(w*ucl/s,n-1)+chi2cdf(w*lcl/s,n-1)-1/t;
    w(floor(t),1)=bisection(1,opw,f,1e-6);
    w(floor(t),2)=bisection(opw,1e6,f,1e-6);
    % 检验结果w对与否
    wj1(floor(t))=1-chi2cdf(w(floor(t),1)*ucl/s,n-1)+chi2cdf(w(floor(t),1)*lcl/s,n-1)-1/t;
    wj2(floor(t))=1-chi2cdf(w(floor(t),2)*ucl/s,n-1)+chi2cdf(w(floor(t),2)*lcl/s,n-1)-1/t;
end
end
function x=bisection(a,b,f,ep)%[a,b]区间，f—函数句柄，ep—最大二分次数
k=1;
eep=b-a;%eep区间长度
while(abs(eep)>ep)
    x=a+eep./2;
    if(f(x).*f(a)<0)
        b=x;
    else
        a=x;
    end
    eep=b-a;
    k=k+1;
end
end

以上是二分法程序