陰解法による数値計算について（1次元拡散方程式）

前提・実現したいこと

1次元拡散方程式を陰解法で解く。

該当のソースコード

python
1import numpy as np
2from matplotlib import pyplot as plt
3from matplotlib import animation as animation
4from pylab import *
5
6
7class Base:
8    def __init__(self, fig, dt, nx, dx):
9        self.dt = float(dt)
10        self.nx = int(nx)
11        self.dx = float(dx)
12        self.dx2 = self.dx * self.dx
13        self.count = 0
14        self.data = np.zeros([self.nx])
15        self.ax = fig.add_subplot(111)
16    def step(self, frame, params):
17        self.count = self.count + 1
18        print(self.count)
19        self.ax.clear()
20        self.ax.grid()
21        data_old = np.copy(self.data)
22        if frame > 0:
23            self.calculation(data_old)
24        self.ax.plot([float(n) * self.dx for n in range(self.nx)], self.data)
25        self.ax.set_xlim(self.xlim)
26        self.ax.set_ylim(self.ylim)
27        self.ax.set_title("t=%.1f" % (float(self.count)*self.dt))
28    def prepare(self):
29        self.prepareData()
30        self.xlim = (0, float(self.nx)*self.dx)
31        self.ylim = (0, np.amax(self.data)*1.2)
32    
33class Sim (Base):
34    def __init__(self, fig, dt, nx, dx, u, a):
35        super().__init__(fig, dt, nx, dx)
36        self.u = float(u)
37        self.a = float(a)
38    def prepareData(self):
39        self.data = np.full(self.nx, 0.0)
40        self.data[4] = 1.0
41        self.data[5] = 1.0
42        self.count = 0
43    def calculation(self, data_old):
44        A = np.zeros((self.nx+1,self.nx+1))
45        F = (self.a * self.dt / self.dx2)
46        A[0,0] = A[10,10] = 2 * F + 1
47        A[0,1] = A[10,9]  = -F
48        for i in range(1, self.nx):
49            if i >= 1 and i < self.nx:
50                A[i,i] = 2 * F + 1
51                A[i,i+1] = A[i,i-1] = -F
52               
53        b = np.zeros(self.nx+1)
54        for ix in range(1,self.nx-1):
55            b[ix] = F * data_old[ix+1] + (1-2*F) * data_old[ix] +F * data_old[ix-1]
56        
57        u = np.zeros(self.nx+1)
58        u[0] =u[10] = 0
59        u =np.linalg.solve(A,b)
60                
61        pass 
62        self.data[0] = self.data[1] 
63        self.data[self.nx-1] = self.data[self.nx-2] 
64
65    
66
67def main():
68    fig = plt.figure()
69    sim = Sim(fig, dt= .1, nx = 10, dx = 1.0, u = 1., a = 1.)
70    ani = animation.FuncAnimation(fig, sim.step,
71            init_func=sim.prepare,
72            frames=range(100),
73            fargs=[sim],
74            blit=False, repeat=False, interval=50)
75    plt.show()
76
77if __name__ == "__main__":
78    main()

試したこと

１次元拡散方程式を陰解法で解こうとしていますが、プログラムでうまく表現できず困っております。
下記ホームページなどを参考にして調べてみたものの、上手く行きませんでした。
https://hplgit.github.io/fdm-book/doc/pub/book/sphinx/._book011.html

初学者のため、基本的なミスや理解の及ばない点があるのは重々承知しておりますが、アドバイス頂けましたら幸いです。

どうぞよろしくお願い致します。