急求二分法修正

玉宇琼儿

程序感觉没错呀，不知为什么，每次运行的时候，b=0,:'( ,坐等解答
Fun[x_] :=
Module[{pw = x},
  sigmaxx = sigmaV;
  sigmayy = sigmaH*Sin[b]^2 + sigmah*Cos[b]^2;
  sigmazz = sigmaH*Cos[b]^2 + sigmah*Sin[b]^2;
  sigmayz = (sigmah - sigmaH)*Cos[b]*Sin[b];
  sigmaxy = 0;
  sigmazx = 0;
  sigmarr = pw;
  pr = 0.25;
  c = 0.95;
  pp = 15*10^6;
  delta = 0.77;
  alpha = 0.55;
  phi = 0.15;
  sigmaH = 45*10^6; sigmah = 35*10^6; sigmaV = 50*10^6; b = 0;
  M = delta*(alpha*(1 - 2*pr)/(1 - pr) - phi)*(pw - pp);
  sigmarr = pw + delta*phi*(pw - pp);
  sigmatt = -pw + (sigmaxx + sigmayy) -
    2*(sigmaxx - sigmayy)*Cos[2*theta] + M;
  sigmaz = -c*pw + sigmazz - pr*(2*(sigmaxx - sigmayy)*Cos[2*theta]) +
     M;
  sigmatz = 2*sigmayz*Cos[b];
  sigmart = 0;
  sigmarz = 0;
  sigma3 =
   0.5*((sigmatt + sigmaz) - Sqrt[(sigmatt - sigmaz)^2 + 4*sigmatz^2]);
  r = FindMinValue[sigma3, {theta}];
  Return[r];
  ]
n = 1;
a0 = 50*10^6;
b0 = 80*10^6;
a = a0;
b = b0;
x = (a + b)/2;
error = 0.00001;
While[Abs[(Fun[x] + 5*10^6)/5/10^6] > error && n < 10, x = (a + b)/2;
If[(Fun[x] + 5*10^6)*(Fun[a] + 5*10^6) < 0, b = x, a = x]; Print[x];
n++]

ggggwhw

我一般看不到别人的代码的,不过我写过牛顿法,二分法,切线法求方程根的代码的