"ZeroDivisionError: float division by zero" error - are my functions declared wrong?

Question

I am trying to write a program which will solve the diffusion equation via the fundamental solution and the representation formula, using the Newton-Cotes composite trapezium to evaluate the integral.

I have attached my work-in-progress code:

import math as m
import matplotlib.pyplot as plt

def initial(xp):
    return (1 / (2 * m.pi * 0.02)**0.5) * (m.exp(-1 * ((xp - 0.2)**2 / (2 * 0.02**2))) + m.exp(-1 * ((xp - 0.8)**2 / (2 * 0.02**2)))) * m.exp(-1 * (1 / (1 - ((2 * xp) - 1)**2)))

def Phi(x,t):
    return (1 / m.sqrt(4 * m.pi * D * t)) * m.exp(- (x**2 / (4 * D * t)))

def integrand(xp,x,t):
    return initial(xp) * Phi((xp-x), t)

def Newton_CT(x,t):
    m = 20
    rest = 0
    for y in range(0, m+1):
        if y == 0:
            donothing = 1
        else:
            rest += integrand((y/m),x,t)
    Qf = 1/(2*m) * (integrand(1,x,t) + 2*rest)
    return Qf

D0 = 10**8
Q = 0.284 * 10**6
T = 1273
R = 8.31

D = D0 * m.exp(-Q / (R * T))

span = 1000
dt = 50
rows = int(span/dt)

time = []
conc = []

for z in range(0, rows+1):
    if z == 0:
        donothing = 1
    else:
        conc.append(Newton_CT(0.5,(z * dt)))
        time.append(z * dt)


fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
plt.grid(True)
ax.plot(time, conc)

The problem comes when I try to run it and I get this string of errors:

  File "...", line 42, in <module>
    conc.append(Newton_CT(0.5,(z * dt)))

  File "...", line 20, in Newton_CT
    rest += integrand((y/m),x,t)

  File "...", line 11, in integrand
    return initial(xp) * Phi((xp-x), t)

  File "...", line 5, in initial
    return (1 / (2 * m.pi * 0.02)**0.5) * (m.exp(-1 * ((xp - 0.2)**2 / (2 * 0.02**2))) + m.exp(-1 * ((xp - 0.8)**2 / (2 * 0.02**2)))) * m.exp(-1 * (1 / (1 - ((2 * xp) - 1)**2)))

ZeroDivisionError: float division by zero

I assume that, somehow, the initial conditions are trying to divide by zero, but I have no idea why that could be happening.

Answer 1

Yes, your functions are declared wrong.

initial is called repeatedly with values for xp increasing in steps of 0.05.

At some point xp is 1, then (1 - ((2 * xp) - 1)**2) will be 0 and the division 1 / (1 - ((2 * xp) - 1)**2) will raise the exception you see.