Use implicit numerical methods to calculate trasient heat temperture.

heat_transient_matrix.py

from mpl_toolkits.mplot3d import axes3d
import numpy as np
import matplotlib.pyplot as plt

n = 15
fo = 0.5

grid = np.zeros(n)
esp = 1e-3
max_iter = None


def init():
    grid[:] = 0
    grid[0] = 0
    grid[-1] = 1


def setMatrix(fo, n):
    a = np.eye(n)
    i = np.arange(n - 2)
    a[i + 1, i] = -fo
    a[i + 1, i + 1] = 1 + 2 * fo
    a[i + 1, i + 2] = -fo
    return getInv(a)


def getInv(a):
    s = a.shape[0]
    ad = np.hstack([a, np.eye(s)])

    for i in range(s):
        # print(i)
        # print(ad)
        ad[i, :] /= ad[i, i]
        j = np.arange(s)
        j = j[j != i]
        ad[j, :] -= ad[i, :] * ad[j, i, np.newaxis]
    return ad[:, s:]


# main
init()
it = 0
history_grid = []
ainv = setMatrix(fo, n)

while True:
    it += 1
    history_grid.append(grid)
    grid = np.dot(ainv, grid)
    if not (np.abs(grid - history_grid[-1]) > esp).any():
        break
    if max_iter and it >= max_iter:
        break

print(it)
print(grid)

ax = plt.gca(projection='3d')
X, Y = np.meshgrid(np.arange(n), np.arange(it))
ax.plot_wireframe(X, Y, np.array(history_grid))

ax.set_xlabel('x')
ax.set_ylabel('timestamp')
ax.set_zlabel('Temperture')
plt.xticks([])
ax.view_init(azim=-140)

plt.savefig("heat_transient_matrix.png")
plt.show()