crank-nicolson格式解决处边值定解问题的python代码

Question

麻烦各位帮我看看代码有啥问题，运行没有结果啊

import numpy as np
h = 1/20
tao = 1/20
i = 1
while 1:
    if abs(np.round(tao,i) - tao) < 10**(-5):
        num_tao = i
        break
    else:
        i += 1
low_space = 0
up_space = 1
low_time = 0
up_time = 1
a = 1
r = a * tao / h**2
x = np.array([
    low_space + h * i for i in range(100)
    if low_space + h * i < 1.0001 * up_space
])
t = np.array([
    low_time + tao * i for i in range(1000)
    if low_time + tao * i < 1.0001 * up_time
])
def exact_value(x,t):
    u = x * (np.e**t + x**2)
    return u
def u_0t(x,t):
    return 0
def u_1t(x,t):
    return 1+np.e**t
def u_x0(x,t):
    return x**3+x
def f_ik(x,t):
    return x * np.e**t - 6 * x
def abcf_computing(x,u_km1,t):
    var = x[1 : -1]
    a = [-r / 2] * len(var)
    b = [1 + r] * len(var)
    c = [-r / 2] * len(var)
    f = []
    for i in range(2,len(var)):
        temp = 0.5 * r * u_km1[i - 1]+(1 - r) * u_km1[i] + 0.5 * r * u_km1[i+1]
        f.append(tao * f_ik(x[i],t - 0.5 * tao)+temp)
        f.insert(0,(1 - r) * u_km1[1] + 0.5 * r * u_km1[2] + tao * f_ik(x[1], - 0.5 * tao + t) + 0.5 * r * (u_0t(0,t) + u_0t(0,t - tao)))
        f.append(u_km1[len(x) - 2] * (1 - r) + u_km1[len(x) - 3] * 0.5 * r + tao * f_ik(x[-2],t - 0.5 * tao) + 0.5 * r *(u_1t(1,t) + u_1t(1,t - tao)))
    return a,b,c,f
def chage_method(x,u_km1,t) :
    a,b,c,f = abcf_computing(x,u_km1,t)
    bbeta = [c[0] / b[0]]
    for i in range(1,len(f) - 1) :
        bbeta.append(c[i] / (b[i] - a[i - 1] * bbeta[i - 1]))
    yy = [f[0] / b[0]]
    for i in range(1,len(f)) :
        yy.append((f[i] - a[i - 1] * yy[i - 1]) / (b[i] - a[i - 1] * bbeta[i - 1]))
    xx = [0 for i in range(len(yy))]
    xx[-1] = yy[-1]
    for i in range(len(yy) - 2,-1,-1) :
        xx[i] = yy[i] - bbeta[i] * xx[i + 1]
    return xx
def crank_nicolson_method() :
    u_i0 = u_x0(x,0)
    result = {'0':u_i0}
    u_ikp1 = list(u_i0.copy())
    for k in range(0,len(t) - 1) :
        u_ik = u_ikp1.copy()
        u_ikp1 = chase_method(x,u_ik,t[k + 1])
        u_ikp1.insert(0,u_0t(0,t[k + 1]))
        u_ikp1.append(u_1t(1,t[k + 1]))
        result[str(np.round((k + 1) * tao,num_tao))] = u_ikp1
    return result
def error() :
    numerical_result = crank_nicolson_method()
    exact_values = {}
    for k in range(len(t)) :
        exact_t = exact_value(x,t[k])
        exact_values[str(np.round(k * tao,num_tao))] = exact_t
    return numerical_result,exact_values
def result_showing() :
    numerical_result,exact_value = error()
    exact_value_04 = []
    result_04 = []
    for i in range(2,11,2) :
        temp1 = numerical_result[str(np.round(i * 0.1,1))]
        temp2 = exact_value[str(np.round(i * 0.1,1))]
        result_04.append(temp1[int(np.round(0.4 / h,0))])
        exact_value_04.append(temp2[int(np.round(0.4 / h,0))])
    errors = abs(np.array(exact_value_04) - np.array(result_04))
    node = [str((0.4,np.round(i * 0.1,1))) for i in range(2,11,2)]
    tplt = "{0:^8}\t{1:^8}\t{2:^8}\t{3:^8}"
    tplt1 = "{0:^10.10}\t{1:^10.8f}\t{2:^10.8f}\t{3:^10.8f}"
    print(tplt.format("节点","数值解","精确解","误差"))
    for i in range(len(exact_value_04)) :
        print(tplt1.format(node[i],result_04[i],exact_value_04[i],errors[i]))
    result_showing()

Answer 1

大致看了 ，主函数没调用，计算得到的结果也没输出