# Power::infy: Infinite expression 1/0. encountered

I am trying to solve the following linear system consisting of four unknowns.

a = 1
b = 1.1
c = 1
w = 100
p = 4500
phi[z_, r_] :=
a5 (8 z^5 - 40 r^2 z^3 + 15 r^4 z) + b5 (2 z^5 - r^2 z^3 - 3 r^4 z) +
a1 (2 z^4 + z^2 r^2 -
r^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)] -
a1 6 z/Sqrt[z^2 + r^2] +
b1 (r^4 + 2 r^2 z^2 +
z^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)]
sr[z_, r_] =
Simplify[D[
v (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {r, 2}], z] -
p w^2 r^2/3 - p w^2 (1 + 2 v) (1 + v) z^2/(6 v (1 - v))]
sz[z_, r_] =
Simplify[D[(2 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], z] +
p w^2 (1 + 3 v) r^2/(6 v)]
trz[z_, r_] =
Simplify[D[(1 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], r]]
NSolve[Integrate[sr[z, a], {z, -c, c}] == 0 &&
Integrate[sr[z, b], {z, -c, c}] == 0 &&
Integrate[trz[c, r] r, {r, a, b}] == 0 &&
Integrate[sz[c, r] r, {r, a, b}] == 0, {a1, b1, a5, b5}]


I remember last week I managed to get a set of numerical answers, but this week Mathematica reports "Power::infy: Infinite expression 1/0. encountered." and keeps running without giving me a result. I don't think there is any infinity point in the function. Can someone tell me how to revise the code to get the answer? Thank you!

• I forget to add that I have defined v to be 0.3. Jun 26, 2019 at 2:35

\$Version

(* "13.1.0 for Mac OS X x86 (64-bit) (June 16, 2022)" *)

Clear["Global*"]


Constants

a = 1;
b = 11/10;
c = 1;
w = 100;
p = 4500;


Functions

phi[z_, r_] =
a5 (8 z^5 - 40 r^2 z^3 + 15 r^4 z) + b5 (2 z^5 - r^2 z^3 - 3 r^4 z) +
a1 (2 z^4 + z^2 r^2 -
r^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)] -
a1 6 z/Sqrt[z^2 + r^2] +
b1 (r^4 + 2 r^2 z^2 +
z^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)] // FullSimplify;

sr[z_, r_] =
FullSimplify[
D[v (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r + D[phi[z, r], {z, 2}]) -
D[phi[z, r], {r, 2}], z] - p w^2 r^2/3 -
p w^2 (1 + 2 v) (1 + v) z^2/(6 v (1 - v))];

sz[z_, r_] =
FullSimplify[
D[(2 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], z] +
p w^2 (1 + 3 v) r^2/(6 v)];

trz[z_, r_] =
FullSimplify[
D[(1 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], r]];


Equations

eqn1 = Integrate[sr[z, a], {z, -c, c}] == 0 // Simplify;

eqn2 = Integrate[sr[z, b], {z, -c, c}] == 0 // Simplify;

eqn3 = Integrate[trz[c, r] r, {r, a, b}] == 0 // Simplify;

eqn4 = Integrate[sz[c, r] r, {r, a, b}] == 0 // Simplify;


The exact solutions are

sol[v_] = Solve[eqn1 && eqn2 && eqn3 && eqn4, {a1, b1, a5, b5}];


The approximate numeric values for v == 0.3 are

sol[0.3]

(* {{a1 -> -6.68076*10^8, b1 -> 8.13029*10^8, a5 -> 4.16484*10^7,
b5 -> 5.18705*10^8}} *)


There are discontinuities for v == 1

sol[1.0] // Quiet

(* {{a1 -> ComplexInfinity, b1 -> Indeterminate,
a5 -> Indeterminate, b5 -> Indeterminate}} *)

Plot[Evaluate[{b1, b5, a5, a1} /. sol[v]], {v, 0, 1.5},
PlotRange -> {-6.8*^8, 8.2*^8},
Frame -> True,
PlotLegends -> Placed[{b1, b5, a5, a1}, {.8, .7}],
Exclusions -> {v == 1},
MaxRecursion -> 5,
GridLines -> {{0.3}, None}]


I made a slight modification to your code and it gives an answer on 12.0.0 + OSX 10.14.5. You may want to make sure that you have nothing defined for any of the internal variables using ClearAll[].

ClearAll[z, v, n, a1, b1, a5, b5];
a = 1;
b = 1.1;
c = 1;
w = 100;
p = 4500;
phi[z_, r_] :=
a5 (8 z^5 - 40 r^2 z^3 + 15 r^4 z) +
b5 (2 z^5 - r^2 z^3 - 3 r^4 z) +
a1 (2 z^4 + z^2 r^2 -
r^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)] -
a1 6 z/Sqrt[z^2 + r^2] +
b1 (r^4 + 2 r^2 z^2 +
z^4) Log[(Sqrt[z^2 + r^2] + z)/(Sqrt[z^2 + r^2] - z)];
sr[z_, r_] =
Simplify[D[
v (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {r, 2}], z] -
p w^2 r^2/3 - p w^2 (1 + 2 v) (1 + v) z^2/(6 v (1 - v))];
sz[z_, r_] =
Simplify[D[(2 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], z] +
p w^2 (1 + 3 v) r^2/(6 v)];
trz[z_, r_] =
Simplify[D[(1 - v) (D[phi[z, r], {r, 2}] + D[phi[z, r], r]/r +
D[phi[z, r], {z, 2}]) - D[phi[z, r], {z, 2}], r]];
NSolve[Integrate[sr[z, a], {z, -c, c}] == 0 &&
Integrate[sr[z, b], {z, -c, c}] == 0 &&
Integrate[trz[c, r] r, {r, a, b}] == 0 &&
Integrate[sz[c, r] r, {r, a, b}] == 0, {a1, b1, a5, b5}]

(* {{a1 -> -((1. (((-267024. +
1.95221*10^-10 I) + (21216.8 -
9.48377*10^-10 I) v + (384367. +
1.45108*10^-9 I) v^2 - (79069.4 +
7.10136*10^-10 I) v^3) (((331.813 +
8.65498*10^-13 I) - (3.41061*10^-13 +
3.38673*10^-13 I) v) (0. + (-10.592 +
11.916 v) (5.13*10^7/(-1. + v) +
5.*10^6/((-1. + v) v) - (
2.63*10^7 v)/(-1. + v))) -
1. (0. + ((2.61056*10^6 - 1.16415*10^-10 I) + (
870187. - 2.91038*10^-11 I)/v) (-10.592 +
11.916 v)) ((-10.592 +
11.916 v) (275.6/(-1. + v) - (
275.6 v)/(-1. + v)) -
52.96 (-(91.12/(-1. + v)) + (149.8 v)/(-1. + v) - (
58.68 v^2)/(-1. + v)))) -
1. (-1. (-1. ((-1.75455 -
1.06581*10^-14 I) - (5.07465 +
7.88861*10^-31 I) v) (1.37707 -
1.67258 v) + ((16.6258 +
0. I) - (13.2092 +
2.84217*10^-14 I) v) (-10.592 +
11.916 v)) ((-10.592 +
11.916 v) (275.6/(-1. + v) - (
275.6 v)/(-1. + v)) -
52.96 (-(91.12/(-1. + v)) + (149.8 v)/(-1. + v) - (
58.68 v^2)/(-1. + v))) + ((331.813 +
8.65498*10^-13 I) - (3.41061*10^-13 +
3.38673*10^-13 I) v) (-1. (1.37707 -
1.67258 v) (-(91.12/(-1. + v)) + (
149.8 v)/(-1. + v) - (
58.68 v^2)/(-1. + v)) + (-10.592 +
11.916 v) (29.6561/(-1. + v) - (
173.856 v)/(-1. + v) + (
144.2 v^2)/(-1. + v)))) (-1. (-52.96 (76. -
36. v) -
200. (-10.592 +
11.916 v)) (0. + ((2.61056*10^6 -
1.16415*10^-10 I) + (870187. - 2.91038*10^-11 I)/
v) (-10.592 + 11.916 v)) + ((331.813 +
8.65498*10^-13 I) - (3.41061*10^-13 +
3.38673*10^-13 I) v) (0. + (-10.592 +
11.916 v) (-3.*10^7 + (
5.*10^6 (1. + v) (1. + 2. v))/((-1. + v) v))))))/((
2.15005*10^8 + 8.27372*10^-7 I)/(-1. +
v) - ((2.30422*10^9 + 7.2686*10^-6 I) v)/(-1. +
v) + ((9.07026*10^9 + 0.0000363434 I) v^2)/(-1. +
v) - ((1.5815*10^10 + 0.0000347323 I) v^3)/(-1. +
v) + ((1.23808*10^10 - 0.0000679472 I) v^4)/(-1. +
v) - ((3.3374*10^9 - 0.000127206 I) v^5)/(-1. +
v) - ((2.0946*10^8 + 0.0000546427 I) v^6)/(-1. +
v) + ((0. + 2.13791*10^-7 I) v^7)/(-1. + v))),
b1 -> -(((4.23724*10^6 -
3.97001*10^6 I) ((-2.68178 +
2.65827 I) + (1.52727*10^16 -
1.63007*10^16 I) v - (3.42115*10^16 -
3.65144*10^16 I) v^2 - (5.0757*10^16 -
5.41736*10^16 I) v^3 + (1.95778*10^17 -
2.08956*10^17 I) v^4 - (1.10133*10^17 -
1.17546*10^17 I) v^5 - (1.58126*10^17 -
1.6877*10^17 I) v^6 + (2.30055*10^17 -
2.4554*10^17 I) v^7 - (1.02746*10^17 -
1.09662*10^17 I) v^8 + (1.5357*10^16 -
1.63907*10^16 I) v^9 - (4.88834*10^14 -
5.21739*10^14 I) v^10 + (1. + 0. I) v^11))/((-1. +
v) v ((3.37708 -
3.27992*10^-14 I) - (0.268331 -
1.44042*10^-14 I) v - (4.86114 - 2.53068*10^-14 I) v^2 +
1. v^3) ((3.87001 -
1.00568*10^15 I) - (33.9987 -
1.07779*10^16 I) v + (169.995 -
4.24259*10^16 I) v^2 - (162.46 -
7.39741*10^16 I) v^3 - (317.821 +
5.79108*10^16 I) v^4 + (595.004 +
1.56106*10^16 I) v^5 - (255.59 -
9.79744*10^14 I) v^6 + (1. + 0. I) v^7))),
a5 -> -(((1.96393*10^-7 -
2.24393*10^-8 I) ((1.7977*10^28 -

1.57338*10^29 I) + (2.1689*10^30 -
1.89826*10^31 I) v - (1.10834*10^31 -
9.70043*10^31 I) v^2 + (2.12511*10^31 -
1.85993*10^32 I) v^3 - (1.52482*10^31 -
1.33455*10^32 I) v^4 - (1.08663*10^31 -
9.51038*10^31 I) v^5 + (3.66366*10^31 -
3.20651*10^32 I) v^6 - (4.01472*10^31 -
3.51376*10^32 I) v^7 + (2.31305*10^31 -
2.02443*10^32 I) v^8 - (6.28849*10^30 -
5.50381*10^31 I) v^9 + (3.8131*10^29 -
3.33729*10^30 I) v^10 + (4.71739*10^28 -
4.12875*10^29 I) v^11 - (9.53322*10^14 -
8.7082*10^14 I) v^12 + (1. + 0. I) v^13))/(((331.813 +
8.65498*10^-13 I) - (3.41061*10^-13 +
3.38673*10^-13 I) v) (-1. +
v) v ((3.37708 -
3.27992*10^-14 I) - (0.268331 -
1.44042*10^-14 I) v - (4.86114 - 2.53068*10^-14 I) v^2 +
1. v^3) ((3.87001 -
1.00568*10^15 I) - (33.9987 -
1.07779*10^16 I) v + (169.995 -
4.24259*10^16 I) v^2 - (162.46 -
7.39741*10^16 I) v^3 - (317.821 +
5.79108*10^16 I) v^4 + (595.004 +
1.56106*10^16 I) v^5 - (255.59 -
9.79744*10^14 I) v^6 + (1. + 0. I) v^7))),
b5 -> ((2.58987*10^-6 -
2.11676*10^-8 I) ((4.39031*10^26 -
5.37159*10^28 I) + (9.38932*10^28 -
1.14879*10^31 I) v - (4.26607*10^29 -
5.21957*10^31 I) v^2 + (4.98212*10^29 -
6.09567*10^31 I) v^3 + (4.92008*10^29 -
6.01977*10^31 I) v^4 - (1.55503*10^30 -
1.9026*10^32 I) v^5 + (9.13331*10^29 -
1.11747*10^32 I) v^6 + (5.72367*10^29 -
7.00296*10^31 I) v^7 - (9.48111*10^29 -
1.16002*10^32 I) v^8 + (4.28735*10^29 -
5.24561*10^31 I) v^9 - (7.30931*10^28 -
8.94301*10^30 I) v^10 + (3.86213*10^27 -
4.72535*10^29 I) v^11 - (9.72616*10^14 -
9.71795*10^14 I) v^12 + (1. + 0. I) v^13))/(((331.813 +
8.65498*10^-13 I) - (3.41061*10^-13 +
3.38673*10^-13 I) v) (-1. + v) v (-0.888889 +
1. v) ((3.37708 -
3.27992*10^-14 I) - (0.268331 -
1.44042*10^-14 I) v - (4.86114 - 2.53068*10^-14 I) v^2 +
1. v^3) ((3.87001 -
1.00568*10^15 I) - (33.9987 -
1.07779*10^16 I) v + (169.995 -
4.24259*10^16 I) v^2 - (162.46 -
7.39741*10^16 I) v^3 - (317.821 +
5.79108*10^16 I) v^4 + (595.004 +
1.56106*10^16 I) v^5 - (255.59 -
9.79744*10^14 I) v^6 + (1. + 0. I) v^7))}} *)

• I forget to add one line of code where I define v to be 0.3. I have also noticed that while Mathematica is able to give me an answer when v is not defined, it reports "Power::infy: Infinite expression 1/0. encountered" when v is defined. Why is this happening? Jun 26, 2019 at 2:45
• Did you click on the stack trace? I ask because I get the same divide by zero problem and the place where it is failing is this: [Integral][Integral]Log[( z + Sqrt[1 + z^2])/(-z + Sqrt[1 + z^2])]/((0. + 1. I) - 1. z)^4 [DifferentialD]z [DifferentialD]z The stack trace may be found by clicking on the button with the 3 dots just to the left of the reported error message. And unfortunately, I don't know why it is failing. Jun 26, 2019 at 5:39
• Problems disappear on my MMA version 8.0, when using infinite precision parameters v = 3/10, b = 11/10  avoiding precision problems. (By the way, use SetDelayed :=  only when necessary, here with phi[z_,r_]:= ` ). Mar 21, 2020 at 21:15