有限差分法求解一维热传导方程

2020-11-01

| 科学计算 | | julia , 有限差分 , 微分方程 , 数值计算 |

# 奇怪：这个库必须放在最前面才能一次加载成功
using Plots
gr()

Plots.GRBackend()

有限差分法求解一维热传导方程

\frac{\partial u}{\partial t}=\lambda \frac{\partial^2 u}{\partial x^2}

具体为(时长T，杆长L)：

\begin{aligned} \frac{\partial u(x,t)}{\partial t}&= \lambda \frac{\partial ^2 u(x,t) }{\partial x^2} \qquad 0\leq t\leq T, 0\leq x \leq L \\ u(x,0)&=4x(L-x)\\ u(0,t)&=0\\ u(L,t)&=0 \end{aligned}

选择合适的步长： $\tau,h$ ，确保 $n=T/\tau,m=L/h$ 都是整数步。

进而有如下递推关系(序号从0开始标注)：

\begin{aligned} u_{i,j+1}&=\frac{\lambda\tau}{h^2}u_{i+1,j}+(1-\frac{2\lambda\tau}{h^2})u_{i,j}+\frac{\lambda\tau}{h^2}u_{i-1,j} \\ & \qquad \qquad \qquad i=1, 2, \dots,(m-1) \qquad j=0,1,\dots,(n-1) \\ u_{i,0}&=4 x_i(L-x_i) \qquad x_i = i h \qquad i=0, 1, \dots,m \\ u_{0,j}&=0 \qquad j=0, 1, \dots,n \\ u_{m,j}&=0 \qquad j=0, 1, \dots,n \end{aligned}

可用Julia实现之（注意代码中数组序号从1开始）：

# 参数确定
λ=1; L=3; T=1; m=30; n=10000;
(h,τ) = (L/m, T/n);

# 网格点坐标序列
t = [(i-1)*τ for i in 1:(n+1)];
x = [(j-1)*h for j in 1:(m+1)];

# 分配空间
u = zeros(m+1,n+1);

# 边界赋值(零值边界不用赋值)
u[:,1] = 4*x.*(L.-x);

# 按递推公式计算
for j = 1:n
    for k = 2:m
        i = m-k+2
        u[i,j+1]=(λ*τ/h^2)*u[i+1,j]+(1-2*λ*τ/h^2)*u[i,j]+(λ*τ/h^2)*u[i-1,j]
    end
end;

plot(t, x, u, xlabel="t", ylabel="x")

《微分方程数值解法（第4版）-李荣华＆刘播-高等教育出版社-2009》