计算线性代数

Question

我尝试求解 1d 不稳定热方程 dT/dt = d^2/dx^2 T(x,t=0) = sin(pi*x) 初始条件和 T(x=0,t)= T(x =1,t) = 0 边界条件，使用 crank-nicolson 方案和 thomas 算法，我正在关闭解决方案，但此图形并不完全满足方程。我找不到我的错在哪里。

`import numpy as np
import matplotlib.pyplot as plt

dt = 0.001
dx = 0.001
n = int(1/dx)
N = int(0.5/dt) # timesteps
R = (2*(dx**2))/dt # parameter for this equation 

# initial zero vectors
x = np.zeros((n+1,1))
b = np.zeros((n+1,N))
a = np.zeros((n+1,N))
c = np.zeros((n+1,N))
d = np.zeros((n+1,N))
F = np.zeros((n+1,N))
Delta = np.zeros((n+1,N))
U = np.zeros((n+1,N))

# points of x axis
for j in range(n+1):
    x[j] = x[0] + (j)*dx

for j in range(N-1):
    for i in range(n+1):
        U[i,0] = np.sin(np.pi*x[i])
        U[n,j] = 0
        U[0,j] = 0

# U = np.round(U,5)

# assign values to p, r, f, b, a, c, d
for j in range(N):
    for i in range(n):
        b[i,j] = 1
        a[i,j] = -(2 + R)
        c[i,j] = 1
for j in range(N):
    for i in range(n):
        d[i,j] = -U[i+1,j] + (2-R)*U[i,j] - U[i-1,j]

# Forward elimination

for j in range(N):
    for i in range(1,n):
        F[0,j] = 0
        Delta[0,j] = 0
        F[i-1,j] = -(b[i-1,j]) / (a[i-1,j] + c[i-1,j]*F[i,j])
        Delta[i-1,j] = (d[i-1,j] - c[i-1,j]*Delta[i,j]) / (a[i-1,j] + c[i-1,j]*F[i,j])

# Back substitution
for j in range(N-1):
    for i in range(n-1,0,-1):
        U[i,j+1] = F[i,j]*U[i-1,j] + Delta[i,j]

plt.plot(x,U)
plt.yticks(np.arange(0,1.1,0.1))
plt.show()`

我得到解决方案图：

Answer 1

我做了一些修改，问题依旧……

在下面的网站上有一个用于 Thomas 算法的函数代码。

“https://folk.ntnu.no/leifh/teaching/tkt4140/._main040.html”

我用过这个功能，我相信它有用。

import numpy as np
import matplotlib.pyplot as plt

dt = 0.05
dx = 0.01
n = int(1/dx) + 1   # number of x positions
N = int(0.5/dt) + 1 # timesteps
R = (2*(dx**2))/dt  # parameter for this equation 

# number of unknown/variable U at each time
# With n points in x, there are (n - 2) variables every time
n_var = n - 2       

# creating/allocating vectors
x = np.linspace(0, 1, n)
b = np.zeros((n_var, N))
a = np.zeros((n_var, N))
c = np.zeros((n_var, N))
d = np.zeros((n_var, N))
F = np.zeros((n_var, N))
Delta = np.zeros((n_var, N))
U = np.zeros((n, N))

# initial condition and boundary conditions
U[1:-1, 0] = np.sin(np.pi*x[1:-1])
U[-1, :] = 0
U[0, :] = 0

# assign values to a, b, c
b[:, :] = 1          # c at tdma function
a[:, :] = -(2 + R)   # b at tdma function
c[:, :] = 1          # a at tdma function


def tdma(a, b, c, d):
    
    """
    Solution of a linear system of algebraic equations with a
        tri-diagonal matrix of coefficients using the Thomas-algorithm.

    Args:
        a(array): an array containing lower diagonal (a[0] is not used)
        b(array): an array containing main diagonal 
        c(array): an array containing lower diagonal (c[-1] is not used)
        d(array): right hand side of the system
    Returns:
        x(array): solution array of the system

    https://folk.ntnu.no/leifh/teaching/tkt4140/._main040.html    
    """
    
    n = len(b)
    x = np.zeros(n)
    
    # elimination:
    
    for k in range(1,n):
        q = a[k]/b[k-1]
        b[k] = b[k] - c[k-1]*q
        d[k] = d[k] - d[k-1]*q
    
    # backsubstitution:
    
    q = d[n-1]/b[n-1]
    x[n-1] = q
    
    for k in range(n-2,-1,-1):
        q = (d[k]-c[k]*q)/b[k]
        x[k] = q
   
    return x


# step by step on time
for j in range(1, N):
    
    # k is equation or variable index (auxiliar indexing)
    # because n_var is different from n
    # i is point/position index (original problem indexing)
    # Use k for a, b, c, d, F and Delta
    # Use i for U
    
    k = np.arange(n_var)
    i = k + 1
    d[k, j] = -U[i+1, j-1] + (2-R)*U[i, j-1] - U[i-1, j-1]
    
    # Its sequence (a, b, c) is different from the function used.
    U[i, j] = tdma(c[k, j], a[k, j], b[k, j], d[k, j])
    

plt.plot(x, U, c='k')
plt.yticks(np.arange(0, 1.1, 0.1))
plt.show()