楼主: 乱云而已
2393 14

[有偿编程] 关于运算精度问题! [推广有奖]

  • 0关注
  • 0粉丝

等待验证会员

已卖:1份资源

高中生

92%

还不是VIP/贵宾

-

威望
0
论坛币
199 个
通用积分
0
学术水平
0 点
热心指数
0 点
信用等级
0 点
经验
231 点
帖子
26
精华
0
在线时间
29 小时
注册时间
2016-4-28
最后登录
2016-5-26

楼主
乱云而已 发表于 2016-5-14 17:04:53 |AI写论文
1800论坛币
下面是我附的代码,给定t值,求w值,一个t对应两个w解,然后将w结果带入方程检验,但是发现有误差。wj2(2315,:)有一个值不为零,想问一下,有谁会用二分法求该方程的解,或者提高该方程的运算精度!

%给定m,n,rou,alpha,opw为p(rou,w)取得最小值,w值,求w两解
clc;clear
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));w=zeros(2315,2);
y=1/(1-chi2cdf(opw*ucl/s,n-1)+chi2cdf(opw*lcl/s,n-1))
for t=1.1:541.1
w(t,:)=fsolve(@(w) 1-chi2cdf(w*ucl/s,n-1)+chi2cdf(w*lcl/s,n-1)-1/t,[1 opw],optimset('Display','off'));
%检验结果w对与否
wj1=zeros(2315,1);
wj2=zeros(2315,1);
wj1(2315,:)=1-chi2cdf(w(t,1)*ucl/s,n-1)+chi2cdf(w(t,1)*lcl/s,n-1)-1/t
wj2(2315,:)=1-chi2cdf(w(t,2)*ucl/s,n-1)+chi2cdf(w(t,2)*lcl/s,n-1)-1/t
end


6W8OV[P)A~6HXJUIKO(2JML.png (7.51 KB)

6W8OV[P)A~6HXJUIKO(2JML.png

关键词:Display fsolve Alpha zeros Optim 二分法

沙发
乱云而已 发表于 2016-5-14 22:08:45
有偿,会的可以谈!

藤椅
cmwei333 发表于 2016-5-15 03:09:18
我写了一个函数,你把 w(t,:)=... 的那一行换成 w(t,:)=f1(t,m,n,rho,alpha,prec); 就可以了,而prec就是你需要的精确度,比如1e-6或者1e-12

f1.zip
下载链接: https://bbs.pinggu.org/a-2032614.html

376 Bytes

需要: 1 个论坛币  [购买]

本附件包括:

  • f1.m

bbs.pinggu.org/forum.php?mod=collection&action=view&ctid=3257
bbs.pinggu.org/forum.php?mod=collection&action=view&ctid=3258
bbs.pinggu.org/forum.php?mod=collection&action=view&ctid=3259

板凳
cmwei333 发表于 2016-5-15 03:11:30
请看我的答复:

https://bbs.pinggu.org/forum.php?mod=redirect&goto=findpost&ptid=4611737&pid=36345645&from^^uid=5975757

报纸
enxizheng 发表于 2016-5-15 11:13:03
貌似你要求解的方程只有一个根, 而不是像你说的一个t对应两个w, 下面是t=1.1时, 对应函数的图像, 可以发现其是单调递减的, 只有一个零点. untitled.jpg

地板
乱云而已 发表于 2016-5-15 15:53:20
enxizheng 发表于 2016-5-15 11:13
貌似你要求解的方程只有一个根, 而不是像你说的一个t对应两个w, 下面是t=1.1时, 对应函数的图像, 可以发现其 ...
函数我证明过了,两个解类似于抛物线形式,还有个解在t比较小时10000以上

7
乱云而已 发表于 2016-5-15 15:55:46
cmwei333 发表于 2016-5-15 03:09
我写了一个函数,你把 w(t,:)=... 的那一行换成 w(t,:)=f1(t,m,n,rho,alpha,prec); 就可以了,而prec就是你 ...
能留下联系方式吗?
不是太明白有些地方!

8
乱云而已 发表于 2016-5-15 15:58:29
enxizheng 发表于 2016-5-15 11:13
貌似你要求解的方程只有一个根, 而不是像你说的一个t对应两个w, 下面是t=1.1时, 对应函数的图像, 可以发现其 ...
两个跟,我证明过,在t比较小时,w的第二个解在1000以上!

9
enxizheng 发表于 2016-5-15 16:33:44
乱云而已 发表于 2016-5-15 15:58
两个跟,我证明过,在t比较小时,w的第二个解在1000以上!
clc;clear
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));w=zeros(2315,2);
y=1/(1-chi2cdf(opw*ucl/s,n-1)+chi2cdf(opw*lcl/s,n-1))

for t=1.1:541.1
    w(floor(t),:)=fsolve(@(w) 1-chi2cdf(w*ucl/s,n-1)+chi2cdf(w*lcl/s,n-1)-1/t,[1 1e5],optimset('Display','off','TolX',1e-15,'MaxIter',1e15,'MaxFunEvals',1e15,'TolFun',1e-15));
    % 检验结果w对与否
    wj1=zeros(2315,1);
    wj2=zeros(2315,1);
    wj1=1-chi2cdf(w(floor(t),1)*ucl/s,n-1)+chi2cdf(w(floor(t),1)*lcl/s,n-1)-1/t
    wj2=1-chi2cdf(w(floor(t),2)*ucl/s,n-1)+chi2cdf(w(floor(t),2)*lcl/s,n-1)-1/t
end

10
enxizheng 发表于 2016-5-15 16:36:55
乱云而已 发表于 2016-5-15 15:58
两个跟,我证明过,在t比较小时,w的第二个解在1000以上!
确实有两个根, 我开始搞错了.
你的算法中第二个初值(opw)给的有点小, 我换成了10000, 结果还挺好的.
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
886.426319511273        20500.6585317666
910.195007401184        19479.8638350741
931.981022318308        18592.4544394032
952.088770334905        17812.0619839549
970.757808373852        17119.0367497982
988.180017947177        16498.4038089188
1004.51145009982        15938.5292493398
1019.88070934880        15430.2226982358
1034.39502283003        14966.1172095453
1048.14472117752        14540.2307685769
1061.20660486635        14147.6498231690
1073.64651280283        13784.2967244464
1085.52130957808        13446.7560356505
1096.88044240670        13132.1428801951
1107.76717480820        12838.0018627932
1118.21957433238        12562.2285083575
1128.27131087319        12303.0075431101
1137.95230750635        12058.7638989533
1147.28927537516        11828.1234984729
1156.30615652488        11609.8814894716
1165.02449305837        11402.9765648995
1173.46373683588        11206.4696692937
1181.64151082901        11019.5265805246
1189.57383099589        10841.4033691420
1197.27529535475        10671.4342520556
1204.75924630293        10509.0213806895
1212.03791030021        10353.6262103604
1219.12251857486        10204.7621699351
1226.02341219182        10061.9884068022
1232.75013326861        9924.90442589658
1239.31150528339        9793.14547591875
1245.71570317296        9666.37856226532
1251.97031526462        9544.29899052862
1258.08239796819        9426.62737144988
1264.05852424836        9313.10691729319
1269.90482670615        9203.50129135117
1275.62703597287        9097.59243229324
1281.23051501540        8995.17884667661
1286.72028986167        8896.07401728293
1292.10107718948        8800.10500598070
1297.37730915024        8707.11121232319
8616.94326781033        8616.94326780494
8529.46204863309        8529.46204862595
8444.53779227202        8444.53779226232
8362.04930521205        8362.04930520650
8281.88325097902        8281.88325097833
8203.93350898463        8203.93350898454
1331.70535507042        8128.10059605520
1336.27705307810        8054.29114414944
1340.77436334844        7982.41742385887
1345.19981617127        7912.39692084499
1349.55581665731        7844.15193955517
1353.84465284763        7777.60925097817
1358.06850317757        7712.69976822658
1362.22944335650        7649.35825090413
1366.32945271665        7587.52303480997
1370.37042007993        7527.13578449710
1374.35414916478        7468.14129766021
1378.28236371760        7410.48714905134
1382.15671196828        7354.12376696860
1385.97877116331        7299.00403537140
1389.75005148778        7245.08319573740
1393.47199982670        7192.31871226805
1397.14600325091        7140.67012413455
1400.77339226577        7090.09891745734
1404.35544384338        7040.56840679187
1407.89338423415        6992.04362529605
1411.38839171410        6944.49122284544
1414.84159886183        6897.87937141550
1418.25409507868        6852.17767713738
1421.62692845484        6807.35709847287
1424.96110886695        6763.38987001367
6720.24943361597        6720.24943145514
6677.91036379152        6677.91036133346
6636.34831794202        6636.34831515512
6595.53997072890        6595.53996757807
6555.46296190473        6555.46295833587
6516.09584563283        6516.09584162145
6477.41804318627        6477.41803867307
6439.40979839189        6439.40979332516
6402.05213597678        6402.05213030869
6365.32682243614        6365.32681610144
6329.21632922470        6329.21632214718
6293.70379814141        6293.70379027589
6258.77300892803        6258.77300017110
6224.40834846354        6224.40833874315
6190.59478205577        6190.59477127803
6157.31782612727        6157.31781424113
6124.56352260908        6124.56350962532
6092.31841506612        6092.31840074093

您需要登录后才可以回帖 登录 | 我要注册

本版微信群
加好友,备注cda
拉您进交流群
GMT+8, 2026-2-3 01:23