Search code examples
pythonscipydifferential-equations

Plotting the solution of diffusion equation for multiple times using SciPy

I want to plot this equation for different times. So time is supposed to be constant, x should vary and then plot y? this equation is the analytic solution of the The Time Dependent Diffusion Equation.

my code so far:

import numpy as np
from scipy.sparse import diags
import scipy as sp
import scipy.sparse
from scipy import special
import matplotlib.pyplot as plt

def Canalytical(intervals, D=1):

    points = 1000

    x=np.linspace(0, 1, intervals+1)
    t=1
    c=np.ones([intervals+1])

    sm = 0
    pos = 0

    for xi in x:

        for i in range(points):

            sm+=sp.special.erfc((1-xi+2*i)/(2*np.sqrt(D*t))) + 
                sp.special.erfc((1+xi+2*i)/(2*np.sqrt(D*t)))

        c[pos] = sm
        pos += 1
        sm = 0

    return c, x

c, xi = Canalytical(intervals=1000) 
plt.plot(xi, c)
plt.show()

I want to plot this equation

Solution

  • The equation in the image is wrong. Plug x = 0 in it, and you'll see it's not zero. The sign in front of the second erfc function should have been -.

    The time t should be passed to Canalytical as a parameter, so the function can be used for multiple values of t.

    Using 1000 terms of the sum is excessive since erfc decays extremely fast at infinity. erfc(10) is about 2e-45, well beyond machine precision, let alone the resolution of the plot.

    Also, consider using vectorization when evaluating functions with NumPy. The entire array x can be passed to the function at once, eliminating the loop. This is what remains: 

    import numpy as np
from scipy import special
import matplotlib.pyplot as plt
def Canalytical(intervals, t=1, D=1):
    points = 1000
    x = np.linspace(0, 1, intervals+1)
    c = np.zeros_like(x)
    for i in range(points):
        c += special.erfc((1-x+2*i)/(2*np.sqrt(D*t))) - special.erfc((1+x+2*i)/(2*np.sqrt(D*t)))
    return x, c
plt.plot(*Canalytical(intervals=1000, t=1)) 
plt.plot(*Canalytical(intervals=1000, t=0.1)) 
plt.plot(*Canalytical(intervals=1000, t=0.01)) 
plt.plot(*Canalytical(intervals=1000, t=0.001)) 
plt.show()

    with the output

    plot