利用向后差分格式计算一维热传导方程初边值问题的数值解的matlab代码

一维热传导方程可以表示为：
$$
\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2}
$$
其中，$u(x,t)$ 为温度场，$\alpha$ 为热扩散系数。

初边值问题可以表示为：
$$
\begin{cases}
\frac{\partial u}{\partial t} = \alpha \frac{\partial^2 u}{\partial x^2} \\
u(x,0) = f(x) \\
u(0,t) = g_1(t) \\
u(L,t) = g_2(t)
\end{cases}
$$
其中，$f(x)$、$g_1(t)$、$g_2(t)$ 分别为初值和边界条件，$L$ 为空间区间长度。

使用向后差分格式离散化时间和空间，可以得到以下公式：
$$
\frac{u_i^{n+1} - u_i^n}{\Delta t} = \alpha \frac{u_{i-1}^{n+1} - 2u_i^{n+1} + u_{i+1}^{n+1}}{\Delta x^2}
$$
其中，$u_i^n$ 表示 $u(x_i, t_n)$ 的近似值，$\Delta t$ 和 $\Delta x$ 分别为时间和空间离散化步长。

将上式整理可得：
$$
u_i^{n+1} = u_i^n + \frac{\alpha \Delta t}{\Delta x^2}(u_{i-1}^{n+1} - 2u_i^{n+1} + u_{i+1}^{n+1})
$$
由于边界条件为 $u(0,t)=g_1(t)$ 和 $u(L,t)=g_2(t)$，因此需要将 $u_0^n$ 和 $u_{N+1}^n$ 固定为 $g_1(t_n)$ 和 $g_2(t_n)$。

以下是 MATLAB 代码实现：

function u = heat_equation_backward(f, g1, g2, alpha, L, N, M)
% f: 初值函数，g1、g2：边界条件
% alpha：热扩散系数，L：空间区间长度
% N：空间分段数，M：时间分段数
dx = L / N; dt = 0.5 * dx^2 / alpha; % 计算 dx 和 dt
x = linspace(0, L, N+1); t = linspace(0, M*dt, M+1);
u = zeros(N+1, M+1); % 存储数值解
% 初始时刻的值
u(:,1) = f(x);
% 边界条件
u(1,:) = g1(t); u(N+1,:) = g2(t);
% 循环计算数值解
for n = 1:M
    for i = 2:N % 注意不用更新边界点
        u(i,n+1) = u(i,n) + alpha*dt/dx^2 * (u(i-1,n+1) - 2*u(i,n+1) + u(i+1,n+1));
    end
end
end

其中，输入参数 f、g1、g2 分别为函数句柄，表示初值函数和边界条件函数。例如，可以定义如下的函数句柄：

f = @(x) sin(pi*x/L); % 初值函数
g1 = @(t) 0; g2 = @(t) 0; % 边界条件为 0

然后调用 `heat_equation_backward` 函数即可计算数值解，例如：

alpha = 0.1; L = 1; N = 100; M = 1000;
u = heat_equation_backward(f, g1, g2, alpha, L, N, M);

该代码计算出了 $[0,L]$ 区间上的温度场 $u(x,t)$ 的数值解。可以使用 `surf` 函数绘制温度场图像：

[X,T] = meshgrid(x,t);
surf(X, T, u');
xlabel('x'); ylabel('t'); zlabel('u');

注意，这里使用了 `meshgrid` 将 $x$ 和 $t$ 向量转化为矩阵，以便于绘制 3D 图像。