一维常系数稳态导热方程的时间显式和时间隐式离散格式的求解程序

一维常系数稳态导热方程的时间显式和时间隐式离散格式如下：

时间显式离散格式：

$\frac{T_i^{n+1} - T_i^n}{\Delta t} = \frac{k}{\rho c_p}\frac{T_{i+1}^n - 2T_i^n + T_{i-1}^n}{(\Delta x)^2}$

其中，$T_i^n$表示时间为$n$时刻位置$i$处的温度，$\Delta t$表示时间步长，$\Delta x$表示空间步长，$k$表示热导率，$\rho$表示密度，$c_p$表示比热容。

时间隐式离散格式：

$\frac{T_i^{n+1} - T_i^n}{\Delta t} = \frac{k}{\rho c_p}\frac{T_{i+1}^{n+1} - 2T_i^{n+1} + T_{i-1}^{n+1}}{(\Delta x)^2}$

我们可以用Matlab编写一个求解程序，如下所示：

时间显式离散格式：

% 常量定义
L = 1; % 杆的长度
nx = 50; % 空间步数
nt = 5000; % 时间步数
dx = L/nx; % 空间步长
dt = 0.001; % 时间步长
k = 1; % 热导率
rho = 1; % 密度
cp = 1; % 比热容

% 初始化温度场
T = zeros(nx+1,1);
T(1) = 0;
T(nx+1) = 1;

% 时间循环
for n = 1:nt
    T_new = zeros(nx+1,1);
    T_new(1) = 0;
    T_new(nx+1) = 1;
    for i = 2:nx
        T_new(i) = T(i) + k*dt/(rho*cp*dx^2)*(T(i+1) - 2*T(i) + T(i-1));
    end
    T = T_new;
end

% 作图
x = 0:dx:L;
plot(x,T);
xlabel('位置');
ylabel('温度');

时间隐式离散格式：

% 常量定义
L = 1; % 杆的长度
nx = 50; % 空间步数
nt = 5000; % 时间步数
dx = L/nx; % 空间步长
dt = 0.001; % 时间步长
k = 1; % 热导率
rho = 1; % 密度
cp = 1; % 比热容

% 初始化温度场
T = zeros(nx+1,1);
T(1) = 0;
T(nx+1) = 1;

% 时间循环
for n = 1:nt
    A = zeros(nx-1,nx-1);
    b = zeros(nx-1,1);
    for i = 1:nx-1
        A(i,i) = 1 + 2*k*dt/(rho*cp*dx^2);
        if i > 1
            A(i,i-1) = -k*dt/(rho*cp*dx^2);
        end
        if i < nx-1
            A(i,i+1) = -k*dt/(rho*cp*dx^2);
        end
        b(i) = T(i+1);
    end
    T_new = [0; A\b; 1];
    T = T_new;
end

% 作图
x = 0:dx:L;
plot(x,T);
xlabel('位置');
ylabel('温度');

这两个程序分别求解了时间显式离散格式和时间隐式离散格式的一维常系数稳态导热方程，并画出了温度随位置的变化图。