差分法求解一维热传导方程matlab程序

teaching matlab

---

$$\begin{aligned} &\frac{\partial u}{\partial t}-a\frac{\partial^2u}{\partial x^2}=0,\\ & u(0,t)=0,\quad u(1,t)=0,\\ & u(x,0)=e^{-\frac{(x-0.25)^2}{0.01}} + 0.1\sin(20\pi x). \end{aligned}$$
其中系数$a=1$. 上述模型的Matlab代码如下:

function pde = model_data()
% MODEL_DATA 模型数据
pde = struct('u_initial',@u_initial,'u_left',@u_left,...
    'u_right',@u_right,'f',@f,'time_grid',@time_grid,...
    'space_grid',@space_grid,'a',@a);
    function [T,tau] = time_grid(NT)
        T = linspace(0,0.1,NT+1);
        tau = 0.1/NT;
    end
    function [X,h] = space_grid(NS)
        X = linspace(0,1,NS+1)';
        h = 1/NS;
    end
    function u = u_initial(x)
        u = exp(-(x-0.025).^2/0.01)+0.1*sin(20*pi*x);
    end
    function u = u_left(t)
        u = zeros(size(t));
    end
    function u = u_right(t)
        u = zeros(size(t));
    end
    function f = f(x,t)
        f = zeros(size(x));
    end
    function a = a()
        a = 1;
    end
end

%% 一维热传导方程有限差分方法主测试脚本 main_test.m
%   依次测试：
%       向前差分
%       向后差分
%       六点对称格式
%   并可视化数值计算结果。
%
% 作者：魏华祎 <weihuayi@xtu.edu.cn> 
pde = model_data(); %模型数据结构体
% 向前差分格式
[X,T,U] = heat_equation_fd1d(100,10000,pde,'forward');
showvarysolution(X,T,U);% 以随时间变化方式显示数值解
showsolution(X,T,U); % 以二元函数方式显示数值解
% 向后差分格式
[X,T,U] = heat_equation_fd1d(100,100,pde,'backward');
showvarysolution(X,T,U);% 以随时间变化方式显示数值解
showsolution(X,T,U); % 以二元函数方式显示数值解
% 六点对称格式，即 Crank-Nicholson 格式
[X,T,U] = heat_equation_fd1d(100,100,pde,'crank-nicholson');
showvarysolution(X,T,U);% 以随时间变化方式显示数值解
showsolution(X,T,U); % 以二元函数方式显示数值解

function [X,T,U] = heat_equation_fd1d(NS,NT,pde,method)
%% HEAT_EQUATION_FD1D 利用有限差分方法计算一维热传导方程
%   
%   输入参数：
%       NS 整型，空间剖分段数
%       NT 整型，时间剖分段数
%       pde 结构体，待求解的微分方程模型的已知数据，
%                  如边界、初始、系数和右端项等条件
%       method 字符串，代表求解所用离散格式
%           F 或 f 或 forward ： 向前差分格式
%           B 或 b 或 backward ： 向后差分格式
%           CN 或 cn 或 crank-nicholson 或 Crank-Nicholson ：
%                      -- 六点对称格式( Crank-Nicholson 格式)
%   输出参数：
%       X 长度为 NS+1 的列向量，空间网格剖分
%       T 长度为 NT+1 的行向量，时间网格剖分
%       U (NS+1)*(NT+1) 矩阵，U(:,i) 表示第 i 个时间层网格部分上的数值解
%
%   作者：魏华祎 <weihuayi@xtu.edu.cn>
[X,h] = pde.space_grid(NS);
[T,tau] = pde.time_grid(NT);
N = length(X);M = length(T);
r = pde.a()*tau/h/h;
if r >= 0.5 && ismember(method,{'F','f','forward'})
    error('时间空间离散不满足向前差分的稳定条件！')
end
U = zeros(N,M);
U(:,1) = pde.u_initial(X);
U(1,:) = pde.u_left(T);
U(end,:) = pde.u_right(T);
switch(method)
    case {'F','f','forward'}
        forward();
    case {'B','b','backward'}
        backward();
    case {'CN','cn','crank-nicholson','Crank-Nicholson'}
        crank_nicholson();
    otherwise
        disp(['Sorry, I do not know your ', method]);
end
%% 向前差分方法
    function forward()
        d = 1 - 2*ones(N-2,1)*r;
        c = ones(N-3,1)*r;
        A = diag(c,-1) + diag(c,1)+diag(d);
        for i = 2:M
            RHS = tau*td.f(X,T(i));
            RHS(2) = RHS(2) + r*U(1,i-1);
            RHS(end-1) = RHS(end-1) + r*U(end,i-1);
            U(2:end-1,i)=A*U(2:end-1,i-1)+ RHS(2:end-1);
        end
    end
%% 向后差分方法
    function backward()
        d = 1 + 2*ones(N-2,1)*r;
        c = -ones(N-3,1)*r;
        A = diag(c,-1) + diag(c,1)+diag(d);    
        for i = 2:M
            RHS = tau*td.f(X,T(i));
            RHS(2) = RHS(2) + r*U(1,i);
            RHS(end-1) = RHS(end-1) + r*U(end,i);
            U(2:end-1,i)=A\(U(2:end-1,i-1)+ RHS(2:end-1));
        end 
    end
%% 六点对称格式， 即 Crank_Nicholson 格式
    function crank_nicholson()
        d1 = 1 + ones(N-2,1)*r;
        d2 = 1 - ones(N-2,1)*r;
        c = 0.5*ones(N-3,1)*r;
        A1 = diag(-c,-1) + diag(-c,1)+diag(d1);  
        A0 = diag(c,-1) + diag(c,1) + diag(d2);
        for i = 2:M
            RHS = tau*td.f(X,T(i));
            RHS(2) = RHS(2) + 0.5*r*(U(1,i)+U(1,i-1));
            RHS(end-1) = RHS(end-1) + ...
                0.5*r*(U(end,i)+U(end,i-1));
            U(2:end-1,i)=A1\(A0*U(2:end-1,i-1)+ RHS(2:end-1));
        end 
    end
end

function showvarysolution(X,T,U)
%%  SHOWVARYSOLUTION  显示数值解随着时间的变化
% 
%  输入参数：
%       X 长度为N的列向量，空间网格剖分
%       T 第度为M的行向量，时间网格剖分
%       U N*M 矩阵，U(:,i) 表示第  i 个时间层网格部分上的数值解
%
%   作者：魏华祎 <weihuayi@xtu.edu.cn>   
M = size(U,2);
figure
xlabel('X');
ylabel('U');
s = [X(1),X(end),min(min(U)),max(max(U))];
axis(s);
for i = 1:M
   plot(X,U(:,i));
   axis(s);
   pause(0.0001);
   title(['T=',num2str(T(i)),' 时刻的温度分布'])
end

function showsolution(X,T,U)
%%  SHOWSOLUTION  以二元函数方式显示数值解
% 
%  输入参数：
%       X 长度为N的列向量，空间网格剖分
%       T 第度为M的行向量，时间网格剖分
%       U N*M 矩阵，U(:,i) 表示第  i 个时间层网格部分上的数值解
%
%   作者：魏华祎 <weihuayi@xtu.edu.cn>   
[x,t] = meshgrid(X,T);
mesh(x,t,U');
xlabel('X');
ylabel('T');
zlabel('U(X,T)');
end

function e = getmaxerror(X,T,U,u_exact)
%%  GETMAXERROR 求最大模误差
%       E(h,\tau) = max_{x_i,t_j}| u_exact(x_i,t_j) - U(i,j)| 
%                 = O( \tau + h^2)
%
% 输入参数：
%       X 长度为 N  的列向量，空间剖分
%       T 长度为 M  的行向量，时间剖分
%       U N*M 的矩阵，U(:,i) 表示第 i 个时间步的数值解
%       u_exact 函数句柄，真解函数
% 输出参数：
%       e 最大模误差
%
% 作者：魏华祎 <weihuayi@xtu.edu.cn>
[x,t] = meshgrid(X,T);
ue = u_exact(x',t');
e = max(max(abs(ue - U)));