-2

I keep getting a divide by zero error in my code and not sure how to fix it.

# coding: utf-8
# In[74]:
import numpy 
import math

Theta_ep = [253.15, 283.15, 313.15]
Theta_e = [[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0],[0.0, 0.0, 0.0]]
p = [1000, 1050, 850, 500, 100]
Lambda = 3.504 
K1 = [0, 0, 0, 0]
K2 = [0, 0, 0, 0]
D = [0, 0, 0, 0,]
R_v = 461          #specific gas constant
Rho = 0.597728
C = 273.15
T_E = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]] 
T_w = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]]
A = 2675.0
e_s = [[0, 0, 0],[0, 0, 0],[0, 0, 0],[0, 0, 0]]
for i in range (0, 4):
    for j in range(0, 3):
        Pi = (p[i+1] / p[0]) ** (1.0 / Lambda)
        K1[i] = -38.5 * (Pi ** (2.0)) + 137.81 * Pi - 53.737
        K2[i] = -4.392 * (Pi ** (2.0)) + 56.831 * Pi - 0.384
        D[i] = (0.1859 * (p[i+1] / p[0]) + 0.6512) ** (-1.0)
        T_E[i][j] = Theta_ep[j] * Pi
        e_s[i][j] = 6.112 * math.exp((17.67*T_E[i][j]) / (T_E[i][j] + 243.5))
        r = 0.622 * (e_s[i][j] / (p[i+1] * 1000.0 - e_s[i][j]))
        Derivative = (4302.65) / ((T_E[i][j] + 243.5) ** 2.0)
        if (C / T_E[i][j]) ** (Lambda) > D[i]:
            T_w[i][j] = -(T_E[i][j] - C - ((A * r) / (1.0 + A * r * Derivative)))            
        elif 1.0 <= (C / T_E[i][j]) ** (Lambda) <= D[i]:
            T_w[i][j] = -(K1[i] - K2[i] * ((C / T_E[i][j]) ** (Lambda)))            
        elif 0.4 <= (C / T_E[i][j]) ** (Lambda) < 1.0:
            T_w[i][j] = ((K1[i] - 1.21) - (K2[i] - 121) * (C / T_E[i][j]) ** (Lambda))            
        elif (C / T_E[i][j]) ** (Lambda) < 0.4:
            T_w[i][j] = ((K1[i] - 2.66) - (K2[i] - 1.21) * (C / T_E[i][j]) ** (Lambda))
        Theta_e[i][j] = T_E[i][j]*((1000.0 / p[i+1]) ** (0.2854 * (1.0-0.28 * 10.0 ** (-3.0 * r)))) * math.exp(((3.376 / T_w[i][j] - 0.00254) * (r * (1.0 + 0.81 * 10.0 ** (-3.0 * r)))))

print ('T_w values:')
print (T_w)

print ('Theta_e values:')
print (Theta_e)

I get a float division by zero error at the first if statement and was wondering if there is a way to fix it?

roganjosh
  • 12,594
  • 4
  • 29
  • 46
  • 4
    Cannot reproduce with the example you've given. Generally we'd ask for a [Minimal, Complete and Verifiable Example](https://stackoverflow.com/help/mcve). This is a lot of code but it does in fact run... you could do with narrowing the issue down to a specific point and it does not throw an error for me in Python 3.6 so it's not verifiable for me. – roganjosh Feb 07 '18 at 20:14
  • It appears to replicate in [Python2.7](https://repl.it/repls/AutomaticEvergreenHammerkop) – TemporalWolf Feb 07 '18 at 20:27
  • @TemporalWolf interesting. 64-bit too? Going to dust off Canopy. – roganjosh Feb 07 '18 at 20:28
  • @roganjosh Dunno, link is to a repl.it which shows the issue, seems to be `T_E[i][j]` is zero at some point in 2.7 when it isn't in 3.6 – TemporalWolf Feb 07 '18 at 20:30
  • @TemporalWolf confirmed in 64-bit 2.7 – roganjosh Feb 07 '18 at 20:31
  • I strongly suspect that the difference between P2.7 and 3.6 is due to integer division somewhere in your code (in P 2.7 you need to cast at least one value to `float`). However, it is taking some effort to debug this for you and since there's so much code I don't wish to keep going. – roganjosh Feb 07 '18 at 20:38
  • 1
    `Pi` in 2.7 evaluates `(850 / 1050 ) ** (1.0 / 3.504)` as `0.0` whereas in 3.6 it is `0.9414772388015634` because `850/1050` is `0` according to integer math. You need to change your `p` values to floats if you want to use 2.7. – TemporalWolf Feb 07 '18 at 20:38

1 Answers1

1

Pi in 2.7 evaluates (850 / 1000 ) ** (1.0 / 3.504) as 0.0 whereas in 3.6 it is 0.9414772388015634 because 850/1000 is 0 according to integer division.

tl;dr: You need to change your p values to floats if you want to use 2.7.

How I figured that out is straight forward:

I wrapped it in a try/except ZeroDivisionError: and then printed variables, tracing the zero back to it's origin:

T_E[i][j] is zero, which is created via T_E[i][j] = Theta_ep[j] * Pi. Theta_ep is never changed and non-zero, so Pi must be zero, which means:

Pi = (p[i+1] / p[0]) ** (1.0 / Lambda) evaluates to zero. Lambda doesn't change and is nonzero, so it must be p[i+1] / p[0] which evaluates to zero.

p[0] == 1000 - In 2.7, integer division rounds down, so anything less than 1000, when divided by 1000, will be zero. The last three values of p all fit that bill.

See repl.it for what I did.

TemporalWolf
  • 7,727
  • 1
  • 30
  • 50