0
$\begingroup$

I am running version 10.1 in Windows 7. The following code gives good results but I have two questions:

  1. Is there a way to speed up the program as I would like to compute more data points?

  2. The code throws a warning or error message

NIntegrate::slwcon: "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. !(*ButtonBox[\">>\", ButtonStyle->\"Link\", ButtonFrame->None, ButtonData:>\"paclet:ref/message/NIntegrate/slwcon\", ButtonNote -> \"NIntegrate::slwcon\"])"

I am looking for suggestions on how to select the options in NIntegrate and FindRoots and also how to set accuracy/precision to make the message goes away.

Mphi[P_, dia_, cover_, fcu1_, fy_, n_, db_, Ab_, fyh_, dsh_, Ash_, 
s_, sg_Integer, tstype_Integer] := 
Module[{Es = 29000, ϵco = 0.002, ϵsp = 0.005, 
Pgoal = 10, Agoal = 10, Maxrecur = 2000, ndata = 10}, 
 ϵy = fy/Es; {ϵsh, ϵsm} = 
If[sg == 1, {14*ϵy, 14*ϵy + 0.14}, 
If[sg == 2, {5*ϵy, 0.12}, Print["error"]]];
rss = ϵsm - ϵsh; fsu = 1.5*fy; 
  m = (fsu/fy)*(((30*rss + 1)^2 - 60*rss - 1)/15/rss^2); 
  fs = 
Piecewise[{{Es*ϵ, \!\(TraditionalForm\`\(-ϵy\) < \
ϵ <= ϵy\)}, 
        {-fy, \!\(TraditionalForm\`\(-ϵsh\) <= \
ϵ < \(-ϵy\)\)}, 
        {(-fy)*((m*(-ϵ - ϵsh) + 
        2)/(60*(-ϵ - ϵsh) + 
        2) + (-ϵ - ϵsh)*
                 ((60 - 
          m)/(2*(30*rss + 
             1)^2))), \!\(TraditionalForm\`\(-ϵsm\) <= \
ϵ < \(-ϵsh\)\)}, 
        {fy, ϵy < ϵ < ϵsh}, {fy*((m*(\
ϵ - ϵsh) + 2)/(60*(ϵ - ϵsh) + 2) + 
               (ϵ - ϵsh)*((60 - 
          m)/(2*(30*rss + 
             1)^2))), \!\(TraditionalForm\`ϵsh <= \
ϵ < ϵsm\)}}]; As = n*Ab; s1 = s - dsh; 
d = dia - cover - db/2; ds = dia - 2*cover + db; 
  rs = dia/2 - cover - db/2; 
tsteel = As/2/Pi/rs; ρs = 4*(Ash/ds/s); ρcc = 
As/((Pi/4)*ds^2); 
ke = If[tstype == 1, (1 - 0.5*s1/ds)/(1 - ρcc), 
If[tstype == 2, (1 - 0.5*s1/ds)^2/(1 - ρcc) , 
 Print["error"]]];
fl = ke*(1/2)*ρs*fyh; Ec = 60.2*Sqrt[fcu1*1000]; 
  Esecu = fcu1/ϵco; ru = Ec/(Ec - Esecu); 
  fcu := 
Piecewise[{{0, ϵ <= 
   0}, {fcu1*(ϵ/ϵco)*(ru/(ru - 
       1 + (ϵ/ϵco)^ru)), 
          \!\(TraditionalForm\`0 < ϵ <= ϵsp\)}, \
{0, ϵ > ϵsp}}]; 
  fcc1 = 
fcu1*(2.254*Sqrt[1 + 7.94*(fl/fcu1)] - 2*(fl/fcu1) - 1.254); 
  ϵcc = ϵco*(1 + 5*(fcc1/fcu1 - 1)); 
Esecc = fcc1/ϵcc; rc = Ec/(Ec - Esecc); 
  ϵcu = ϵsp + 
1.4*ρs*(fyh/fcc1)*ϵsm; 
  fcc := 
Piecewise[{{0, ϵ <= 
   0}, {fcc1*(ϵ/ϵcc)*(rc/(rc - 
       1 + (ϵ/ϵcc)^rc)), 
          \!\(TraditionalForm\`0 < ϵ <= ϵcu\)}, \
{0, ϵ > ϵcu}}]; 
  Acover = 
ImplicitRegion[(ds/2)^2 <= x^2 + y^2 <= (dia/2)^2, {x, y}]; 
  Acore = ImplicitRegion[x^2 + y^2 <= (ds/2)^2, {x, y}]; 
  Asteel = ImplicitRegion[x^2 + y^2 == rs^2, {x, y}]; 
 modratio = Es/Ec; Ag = (Pi/4)*dia^2; 
  ϵaxial = 
 P/Ec/Ag; Δϵ = (ϵcu - \
ϵaxial)/ndata; 
corestrain = 
Range[ϵaxial + Δϵ, ϵcu, \
Δϵ]; 
  iterate[ϵbstart_, ϵcore_] := 
Module[{}, Ccover[ϵb_?NumericQ] := 

 NIntegrate[
  fcu /. ϵ -> ϵcore - ((ϵcore - \
 ϵb)/(d - cover + 0.5*dsh))*ycore /. ycore -> ds/2 - y, 
           \!\(TraditionalForm\`{x, y} ∈ Acover\)]; 
 Ccore[ϵb_?NumericQ] := 

 NIntegrate[
  fcc /. ϵ -> ϵcore - ((ϵcore - \
 ϵb)/(d - cover + 0.5*dsh))*ycore /. ycore -> ds/2 - y, 
          \!\(TraditionalForm\`{x, y} ∈ Acore\)]; 
 Cs[ϵb_?NumericQ] := 

 NIntegrate[(fs - fcc)*
     tsteel /. ϵ -> ϵcore - ((ϵcore - \
ϵb)/(d - cover + 0.5*dsh))*ycore /. ycore -> ds/2 - y, 
           \!\(TraditionalForm\`{x, y} ∈ Asteel\), 
  MaxRecursion -> Maxrecur, PrecisionGoal -> Pgoal, 
           AccuracyGoal -> Agoal]; ϵb /. 
 FindRoot[
  Ccover[ϵb] + Ccore[ϵb] + Cs[ϵb] == P, 
           {ϵb, ϵbstart}, 
  MaxIterations -> 100]]; 
 soln = FoldList[iterate[#1, #2] & , ϵaxial, corestrain]; 
  strains = Transpose[{corestrain, Rest[soln]}]; 
  ϕ = ((strains[[#1, 1]] - strains[[#1, 2]])/(d - cover + 
     0.5*dsh) & ) /@ Range[Length[corestrain]];     
 munconfined = (NIntegrate[
   fcu*y /. ϵ -> #1[[
        1]] - ((#1[[1]] - #1[[2]])/(d - cover + 0.5*dsh))*ycore /. 
    ycore -> ds/2 - y, \!\(TraditionalForm\`{x, y} ∈ 
      Acover\), MaxRecursion -> Maxrecur, 
           PrecisionGoal -> Pgoal, AccuracyGoal -> Agoal] & ) /@ 
  strains;     
mconfined = (NIntegrate[
   fcc*y /. ϵ -> #1[[
        1]] - ((#1[[1]] - #1[[2]])/(d - cover + 0.5*dsh))*ycore /. 
    ycore -> 
     ds/2 - y, \!\(TraditionalForm\`{x, y} ∈ Acore\), 
   MaxRecursion -> Maxrecur, 
           PrecisionGoal -> Pgoal, AccuracyGoal -> Agoal] & ) /@ 
strains;     
msteel = (NIntegrate[(fs - fcc)*tsteel*
      y /. ϵ -> #1[[
        1]] - ((#1[[1]] - #1[[2]])/(d - cover + 0.5*dsh))*ycore /. 
    ycore -> ds/2 - y, \!\(TraditionalForm\`{x, y} ∈ 
      Asteel\), MaxRecursion -> Maxrecur, 
           PrecisionGoal -> Pgoal, AccuracyGoal -> Agoal] & ) /@ 
strains; 
  moment = munconfined + mconfined + msteel; 
neutralaxis = corestrain[[#1]]/ϕ[[#1]] + cover - 0.5*dsh & /@
      Range@Length@corestrain; 
data = {ϕ, moment, corestrain, Rest[soln], 
 neutralaxis}\[Transpose]; data]

To run the code:

Mphi[350, 24, 1, 5, 48.5, 10, 1, 0.79, 45, 0.375, 0.11, 3.5, 1, 2]
$\endgroup$
  • $\begingroup$ Which integral exactly issues the warning? $\endgroup$ – Michael E2 Jun 5 '15 at 21:59
  • $\begingroup$ I am not sure which integral gave the warning. I am hoping you might be able to help me find out :) $\endgroup$ – user11946 Jun 6 '15 at 0:13
  • 1
    $\begingroup$ I can help: use Print to mark different points in your code. $\endgroup$ – Gregory Rut Jun 6 '15 at 8:37
  • $\begingroup$ The variables ϵ does not have a numeric value, and that makes everything turn in to a huge mess. $\endgroup$ – Ted Ersek Feb 2 '17 at 1:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.