能不能别白嫖啊٩(๑❛ᴗ❛๑)۶
目录
一、前言
1.本人在做建模比赛2020年A题时遇到热传导方程，然后我学习并总结了用有限差分法求解此类偏微分方程（抛物型方程）的解法。
2.本文借鉴《偏微分方程数值解法——孙志忠》，这本书写得浅显易懂。
3.本文不讲复杂的推理过程，只讲方法与结果。
4.相较其他博客而言（绝大多数都是用matlab），我的特点使用Python实现程序的
5.如有问题可在评论区讨论或是私信我
二、问题与定义
1.问题
考虑一维非齐次热传导方程的定解问题
2.定义
对x轴[0,1]区间作m等分，将t轴[0,T]区间作n等分，并记h=1/m，τ=T/n，分别称h和τ为空间步长和时间步长。
记
为网格比。
记 
网格剖分如下所示
定义网格函数
其中
是的近似
三、求解格式
1、向前欧拉格式（古典显格式）
（1）差分格式
可以看出k+1层结点可由第k层得出，因此可以简单循环即可
（2）收敛性与稳定性
当r<=1/2时是收敛与稳定的，但是当r>1/2时是不稳定的，也不一定收敛。
2.向后欧拉格式（古典隐格式）
（1）差分格式
可以看出第k层的结点由k+1层结点得出，不能直接迭代求解，必须解方程组
将如上差分格式写为矩阵形式（注意不写的地方为0）
只需要在每一个时间层解这个三对角线性方程组即可，用AX=B代替以上方程组以方便叙述
有两种方法：第一种是直接将方程组两端乘上A逆即可得到X；第二种是用三对角线性方程组的特殊求解方法：追赶法来求解，当矩阵很大的时候速度比求逆矩阵要快得多，方法并不难，我放在最后
（2）收敛性与稳定性
对于任意步长比r必收敛与稳定
3.Crank-Nicolson格式
（1）差分格式
同样是在每一个时间层解一个三对角线性方程组即可
（2）收敛性与稳定性
对于任意步长比r必收敛与稳定
四、追赶法求解三对角线性方程组
追赶法
五、实例+代码
以下题目的结果均保存至excel中方便观察
1.题目1（《偏微分方程数值解法——孙志忠》书中题目）
时间步长与空间步长均取0.1
代码：
import numpy as np import pandas as pd import datetime start_time = datetime.datetime.now() np.set_printoptions(suppress=True) def left_boundary(t): # 左边值 return np.exp(t) def right_boundary(t): # 右边值 return np.exp(t + 1) def initial_T(x_max, t_max, delta_x, delta_t, m, n): # 给温度T初始化 T = np.zeros((n + 1, m + 1)) for i in range(m + 1): # 初值 T[0, i] = np.exp(i * delta_x) for i in range(1, n + 1): # 注意不包括T[0,0]与T[0,-1] T[i, 0] = left_boundary(i * delta_t) # 左边值 T[i, -1] = right_boundary(i * delta_t) # 右边值 return T # 一、古典显格式 def one_dimensional_heat_conduction1(T, m, n, r): # 可以发现当r>=0.5时就发散了 for k in range(1, n + 1): # 时间层 for i in range(1, m): # 空间层 T[k, i] = (1 - 2 * r) * T[k - 1, i] + r * (T[k - 1, i - 1] + T[k - 1, i + 1]) return T.round(6) # 二、古典隐格式（乘逆矩阵法） def one_dimensional_heat_conduction2(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + 2 * r) + np.eye(m - 1, k=1) * (-r) + np.eye(m - 1, k=-1) * (-r) a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): T[k, 1:-1] = np.linalg.inv(A) @ F # 左乘A逆 return T.round(6) # 三、古典隐格式（追赶法） def one_dimensional_heat_conduction3(T, m, n, r): a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k, 1:-1] = x return T.round(6) # 四、Crank-Nicolson（乘逆矩阵法） def one_dimensional_heat_conduction4(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + r) + np.eye(m - 1, k=1) * (-r * 0.5) + np.eye(m - 1, k=-1) * (-r * 0.5) C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r / 2 * (T[k, 0] + T[k + 1, 0]) F[-1] = r / 2 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F T[k + 1, 1:-1] = np.linalg.inv(A) @ F return T.round(6) # 五、Crank-Nicolson（追赶法） def one_dimensional_heat_conduction5(T, m, n, r): C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) a = np.ones(m - 1) * (-0.5 * r) a[0] = 0 b = np.ones(m - 1) * (1 + r) c = np.ones(m - 1) * (-0.5 * r) c[-1] = 0 for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r * 0.5 * (T[k, 0] + T[k + 1, 0]) F[-1] = r * 0.5 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k + 1, 1:-1] = x return T.round(6) def exact_solution(T, m, n, r, delta_x, delta_t): # 偏微分方程精确解 for i in range(n + 1): for j in range(m + 1): T[i, j] = np.exp(i * delta_t + j * delta_x) return T.round(6) a = 1 # 热传导系数 x_max = 1 t_max = 1 delta_x = 0.1 # 空间步长 delta_t = 0.1 # 时间步长 m = int((x_max / delta_x).__round__(4)) # 长度等分成m份 n = int((t_max / delta_t).__round__(4)) # 时间等分成n份 t_grid = np.arange(0, t_max + delta_t, delta_t) # 时间网格 x_grid = np.arange(0, x_max + delta_x, delta_x) # 位置网格 r = (a * delta_t / (delta_x 2)).__round__(6) # 网格比 T = initial_T(x_max, t_max, delta_x, delta_t, m, n) print('长度等分成{}份'.format(m)) print('时间等分成{}份'.format(n)) print('网格比=', r) p = pd.ExcelWriter('有限差分法-一维热传导-题目1.xlsx') T1 = one_dimensional_heat_conduction1(T, m, n, r) T1 = pd.DataFrame(T1, columns=x_grid, index=t_grid) # colums是列号，index是行号 T1.to_excel(p, '古典显格式') T2 = one_dimensional_heat_conduction2(T, m, n, r) T2 = pd.DataFrame(T2, columns=x_grid, index=t_grid) # colums是列号，index是行号 T2.to_excel(p, '古典隐格式（乘逆矩阵法）') T3 = one_dimensional_heat_conduction3(T, m, n, r) T3 = pd.DataFrame(T3, columns=x_grid, index=t_grid) # colums是列号，index是行号 T3.to_excel(p, '古典隐格式（追赶法）') T4 = one_dimensional_heat_conduction4(T, m, n, r) T4 = pd.DataFrame(T4, columns=x_grid, index=t_grid) # colums是列号，index是行号 T4.to_excel(p, 'Crank-Nicolson格式（乘逆矩阵法）') T5 = one_dimensional_heat_conduction5(T, m, n, r) T5 = pd.DataFrame(T5, columns=x_grid, index=t_grid) # colums是列号，index是行号 T5.to_excel(p, 'Crank-Nicolson格式（追赶法）') T6 = exact_solution(T, m, n, r, delta_x, delta_t) T6 = pd.DataFrame(T6, columns=x_grid, index=t_grid) # colums是列号，index是行号 T6.to_excel(p, '偏微分方程精确解') p.save() end_time = datetime.datetime.now() print('运行时间为', (end_time - start_time)) 书中Crank-Nicolson的数值解
可以看出与我的程序答案完全一致。
其余格式的结果与我的程序答案都是一样的。
结果分析大家就自己做了，我就不赘述了。
2.题目2
（1）理论分析
由生活实际知，温度不会无限制增大，当时间趋于一定的值时，温度将会稳定。
而由一维热传导方程知，
温度与距离将会呈线性关系，若以杆最左端为x=0则有：y=175+4x(x=5时y=195)
在这个题目中，f(x,t)=0，差分格式用矩阵比较容易表达与求解
如果f(x,t)≠0，那么可以自己在程序中修改矩阵的值
（2）代码
import numpy as np import pandas as pd import datetime start_time = datetime.datetime.now() np.set_printoptions(suppress=True) def left_boundary(t): # 左边值 T = 25 + 15 * t if T < 175: return T else: return 175 def right_boundary(t): # 右边值 T = 25 + 17 * t if T < 195: return T else: return 195 def initial_T(x_max, t_max, delta_x, delta_t, m, n): # 给温度T初始化 T = np.zeros((n + 1, m + 1)) for i in range(m + 1): # 初值 T[0, i] = 25 for i in range(1, n + 1): # 注意不包括T[0,0]与T[0,-1] T[i, 0] = left_boundary(i * delta_t) # 左边值 T[i, -1] = right_boundary(i * delta_t) # 右边值 return T # 一、古典显格式 def one_dimensional_heat_conduction1(T, m, n, r): # 可以发现当r>=0.5时就发散了 for k in range(1, n + 1): # 时间层 for i in range(1, m): # 空间层 T[k, i] = (1 - 2 * r) * T[k - 1, i] + r * (T[k - 1, i - 1] + T[k - 1, i + 1]) return T.round(6) # 二、古典隐格式（乘逆矩阵法） def one_dimensional_heat_conduction2(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + 2 * r) + np.eye(m - 1, k=1) * (-r) + np.eye(m - 1, k=-1) * (-r) a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): T[k, 1:-1] = np.linalg.inv(A) @ F # 左乘A逆 return T.round(6) # 三、古典隐格式（追赶法） def one_dimensional_heat_conduction3(T, m, n, r): a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k, 1:-1] = x return T.round(6) # 四、Crank-Nicolson（乘逆矩阵法） def one_dimensional_heat_conduction4(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + r) + np.eye(m - 1, k=1) * (-r * 0.5) + np.eye(m - 1, k=-1) * (-r * 0.5) C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r / 2 * (T[k, 0] + T[k + 1, 0]) F[-1] = r / 2 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F T[k + 1, 1:-1] = np.linalg.inv(A) @ F return T.round(6) # 五、Crank-Nicolson（追赶法） def one_dimensional_heat_conduction5(T, m, n, r): C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) a = np.ones(m - 1) * (-0.5 * r) a[0] = 0 b = np.ones(m - 1) * (1 + r) c = np.ones(m - 1) * (-0.5 * r) c[-1] = 0 for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r * 0.5 * (T[k, 0] + T[k + 1, 0]) F[-1] = r * 0.5 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k + 1, 1:-1] = x return T.round(6) a = 0.5 # 热传导系数 x_max = 5 t_max = 100 delta_x = 0.1 # 空间步长 delta_t = 0.1 # 时间步长 m = int((x_max / delta_x).__round__(4)) # 长度等分成m份 n = int((t_max / delta_t).__round__(4)) # 时间等分成n份 t_grid = np.arange(0, t_max + delta_t, delta_t) # 时间网格 x_grid = np.arange(0, x_max + delta_x, delta_x) # 位置网格 r = (a * delta_t / (delta_x 2)).__round__(6) # 网格比 T = initial_T(x_max, t_max, delta_x, delta_t, m, n) print('长度等分成{}份'.format(m)) print('时间等分成{}份'.format(n)) print('网格比=', r) p = pd.ExcelWriter('有限差分法-一维热传导-题目2.xlsx') T1 = one_dimensional_heat_conduction1(T, m, n, r) T1 = pd.DataFrame(T1, columns=x_grid, index=t_grid) # colums是列号，index是行号 T1.to_excel(p, '古典显格式') T2 = one_dimensional_heat_conduction2(T, m, n, r) T2 = pd.DataFrame(T2, columns=x_grid, index=t_grid) # colums是列号，index是行号 T2.to_excel(p, '古典隐格式（乘逆矩阵法）') T3 = one_dimensional_heat_conduction3(T, m, n, r) T3 = pd.DataFrame(T3, columns=x_grid, index=t_grid) # colums是列号，index是行号 T3.to_excel(p, '古典隐格式（追赶法）') T4 = one_dimensional_heat_conduction4(T, m, n, r) T4 = pd.DataFrame(T4, columns=x_grid, index=t_grid) # colums是列号，index是行号 T4.to_excel(p, 'Crank-Nicolson格式（乘逆矩阵法）') T5 = one_dimensional_heat_conduction5(T, m, n, r) T5 = pd.DataFrame(T5, columns=x_grid, index=t_grid) # colums是列号，index是行号 T5.to_excel(p, 'Crank-Nicolson格式（追赶法）') p.save() end_time = datetime.datetime.now() print('运行时间为', (end_time - start_time)) （3）运行结果与分析
可以发现显格式已经不收敛了，其余格式都收敛；
最终的温度都稳定为y=175+4x，符合理论分析。
3.特殊的例子：题目3
这里仅仅是修改了边值条件，但是情况却大有不同了！
（1）代码
注：结果会存到“有限差分法-一维热传导.xlsx”中，方便大家观察
import numpy as np import pandas as pd import datetime start_time = datetime.datetime.now() np.set_printoptions(suppress=True) def left_boundary(t): # 左边 return 175 def right_boundary(t): # 右边值 return 195 def initial_T(x_max, t_max, delta_x, delta_t, m, n): # 给温度T初始化 T = np.zeros((n + 1, m + 1)) for i in range(m + 1): # 初值 T[0, i] = 25 for i in range(1, n + 1): # 注意不包括T[0,0]与T[0,-1] T[i, 0] = left_boundary(i * delta_t) # 左边值 T[i, -1] = right_boundary(i * delta_t) # 右边值 return T # 一、古典显格式 def one_dimensional_heat_conduction1(T, m, n, r): # 可以发现当r>=0.5时就发散了 for k in range(1, n + 1): # 时间层 for i in range(1, m): # 空间层 T[k, i] = (1 - 2 * r) * T[k - 1, i] + r * (T[k - 1, i - 1] + T[k - 1, i + 1]) return T.round(6) # 二、古典隐格式（乘逆矩阵法） def one_dimensional_heat_conduction2(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + 2 * r) + np.eye(m - 1, k=1) * (-r) + np.eye(m - 1, k=-1) * (-r) a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): T[k, 1:-1] = np.linalg.inv(A) @ F # 左乘A逆 return T.round(6) # 三、古典隐格式（追赶法） def one_dimensional_heat_conduction3(T, m, n, r): a = np.ones(m - 1) * (-r) a[0] = 0 b = np.ones(m - 1) * (1 + 2 * r) c = np.ones(m - 1) * (-r) c[-1] = 0 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) for k in range(1, n + 1): # 时间层range(1, n + 1) F[0] = T[k - 1, 1] + r * T[k, 0] F[-1] = T[k - 1, m - 1] + r * T[k, m] y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 2): # 空间层 F[i] = T[k - 1, i + 1] # 给F赋值 for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k, 1:-1] = x return T.round(6) # 四、Crank-Nicolson（乘逆矩阵法） def one_dimensional_heat_conduction4(T, m, n, r): A = np.eye(m - 1, k=0) * (1 + r) + np.eye(m - 1, k=1) * (-r * 0.5) + np.eye(m - 1, k=-1) * (-r * 0.5) C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r / 2 * (T[k, 0] + T[k + 1, 0]) F[-1] = r / 2 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F T[k + 1, 1:-1] = np.linalg.inv(A) @ F return T.round(6) # 五、Crank-Nicolson（追赶法） def one_dimensional_heat_conduction5(T, m, n, r): C = np.eye(m - 1, k=0) * (1 - r) + np.eye(m - 1, k=1) * (0.5 * r) + np.eye(m - 1, k=-1) * (0.5 * r) a = np.ones(m - 1) * (-0.5 * r) a[0] = 0 b = np.ones(m - 1) * (1 + r) c = np.ones(m - 1) * (-0.5 * r) c[-1] = 0 for k in range(0, n): # 时间层 F = np.zeros(m - 1) # m-1个元素，索引0~(m-2) F[0] = r * 0.5 * (T[k, 0] + T[k + 1, 0]) F[-1] = r * 0.5 * (T[k, m] + T[k + 1, m]) F = C @ T[k, 1:m] + F y = np.zeros(m - 1) beta = np.zeros(m - 1) x = np.zeros(m - 1) y[0] = F[0] / b[0] d = b[0] for i in range(1, m - 1): beta[i - 1] = c[i - 1] / d d = b[i] - a[i] * beta[i - 1] y[i] = (F[i] - a[i] * y[i - 1]) / d x[-1] = y[-1] for i in range(m - 3, -1, -1): x[i] = y[i] - beta[i] * x[i + 1] T[k + 1, 1:-1] = x return T.round(6) a = 0.5 # 热传导系数 x_max = 5 t_max = 100 delta_x = 0.1 # 空间步长 delta_t = 0.1 # 时间步长 m = int((x_max / delta_x).__round__(4)) # 长度等分成m份 n = int((t_max / delta_t).__round__(4)) # 时间等分成n份 t_grid = np.arange(0, t_max + delta_t, delta_t) # 时间网格 x_grid = np.arange(0, x_max + delta_x, delta_x) # 位置网格 r = (a * delta_t / (delta_x 2)).__round__(6) # 网格比 T = initial_T(x_max, t_max, delta_x, delta_t, m, n) print('长度等分成{}份'.format(m)) print('时间等分成{}份'.format(n)) print('网格比=', r) p = pd.ExcelWriter('有限差分法-一维热传导-题目3.xlsx') T1 = one_dimensional_heat_conduction1(T, m, n, r) T1 = pd.DataFrame(T1, columns=x_grid, index=t_grid) # colums是列号，index是行号 T1.to_excel(p, '古典显格式') T2 = one_dimensional_heat_conduction2(T, m, n, r) T2 = pd.DataFrame(T2, columns=x_grid, index=t_grid) # colums是列号，index是行号 T2.to_excel(p, '古典隐格式（乘逆矩阵法）') T3 = one_dimensional_heat_conduction3(T, m, n, r) T3 = pd.DataFrame(T3, columns=x_grid, index=t_grid) # colums是列号，index是行号 T3.to_excel(p, '古典隐格式（追赶法）') T4 = one_dimensional_heat_conduction4(T, m, n, r) T4 = pd.DataFrame(T4, columns=x_grid, index=t_grid) # colums是列号，index是行号 T4.to_excel(p, 'Crank-Nicolson格式（乘逆矩阵法）') T5 = one_dimensional_heat_conduction5(T, m, n, r) T5 = pd.DataFrame(T5, columns=x_grid, index=t_grid) # colums是列号，index是行号 T5.to_excel(p, 'Crank-Nicolson格式（追赶法）') p.save() end_time = datetime.datetime.now() print('运行时间为', (end_time - start_time)) （2）运行结果与分析
显格式仍不收敛。
注意！Crank-Nicolson格式虽然最终收敛于y=175+4x，但是可以发现在最初的1秒内出现了数值伪振荡，这是因为我们给定的边值不合理导致的。通俗的讲，试想在生活实际中那能让一根杆在极短时间从25摄氏度跃变至175摄氏度，这合理吗？并不合理，这么大的变化率导致CN格式出现了数值伪振荡。并且只有当网格比稍大的时候才会出现这种情况，网格比小时仍然可用CN格式
版权声明：本文内容由互联网用户自发贡献，该文观点仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请联系我们举报，一经查实，本站将立刻删除。
发布者：全栈程序员-站长，转载请注明出处：https://javaforall.net/227449.html原文链接：https://javaforall.net
