偏微分方程实验报告二
实验题目：利用向前差分求解如下一维热传导方程，并观察最大规模误差变化情况：
?u?2u?a2?f(x,t),?t?xu(0,t)?uL(x),u(1,t)?uR(x), u(x,0)?u0(x),x?[0,1],t?[0,0.1].其中a?1;真解为u(x,t)?sin(2?x)e10t 最大模误差定义如下：
E?maxu(xi,tj)?U(i,j)
xi,tj
实现算法：对于初边值问题（抛物方程）
?u?2u?a2?f(x),(x,t)?G,?t?x (1) u(x,0)??(x),0?x?l,u(0,t)?u(l,t)?0,0?t?T,其相应的有限差分法的算法如下：
1.对求解区域做网格剖分，得到计算网格.取空间步长和时间步长为h?1T,??，分别对NM空间变量X所属的区间[0,l]和时间变量t所属的区间[0,T]做如下均匀剖分：
0?x0?x1?...?xN?l,0?t0?t1?...tM?T
其中xj?jh,tk?k?.
用两族平行直线x?xj(0,1,...,N)和t?tk(k?0,1,...,M)将矩形域G分割成矩形网格. 2.对微分方程中的各阶导数进行差分离散，得到差分方程.用uj表示差分解在网格节点处的分量，下面采用逐层计算的思想建立差分格式. （xj,tk）设从第0个时间层到第k?0个时间层的差分解分量
kuij,i?0,...,k;j?0,...,N
已经求得.
下面建立第k?1个时间层（简称当前层）上的向前差分方程，对当前层上的任意内节点
(xj,tk?1)，规定其差分格式涉及的模板点为
对处的微分方程 （xj,tk）?u?2u?a2?t(xj,tk)?x(x中的导数用差商近似：
?f(xj) （1）
j,tk)?u(xj,tk)?t?2u(xj,tk)?x2?u(xj,tk?1)?u(xj,tk)（向前差商）??u(xj?1,tk)?2u(xj,tk)?u(xj?1,tk)h2（2）
(二阶中心差商)记fj?f(xj)，将（2）代入（1），则有关于的向前差分格式 （xj,tk?1）?1uk?ukjj?在（3）两边同乘?，整理可得
?akkukj?1?2uj?uj?1h2?fj （3）
?1kkuk?rukjj?1?(1?2r)uj?ruj?1??fj （4）
其中，r?a?h2称为网格比（简称网比）.
3.根据初边值条件，进行初边值处理.由(4)可得
?1kkuk?ruk;k?0,...,M?1jj?1?(1?2r)uj?ruj?1??fj,j?1,...,N?1u0j?u(xj,0)??(xj),j?1,...,N?1kku0?u(0,tk)?0,uN?u(l,tk)?0,k?0,1,...,M（5）
令N?1维（第k个网格层）差分解向量和右端向量
kkTTUk?(u1k,u2,...,uN?1),F?(f1,f2,...,fN?1)
则向前差分格式（4）的矩阵表示为：
Uk?1?A0Uk??F
其中
4.最后求解线性代数方程组，得到数值解向量Un?1
.
程序代码：
function [X,T,U]=heat_equation_fd1d(NS,NT,pde) %% HEAT_EQUATION_FD1D 利用有限差分方法计算一维热传导方程 % 输入参数：
% NS _x0012_整型，空间剖分段数 % NT _x0012_整型，时间剖分段数
% pde 结构体，待求解的微分方程模型的已知数据，如边界、初始、系数和右端项等条件 % 输出参数
% X 长度为 NS +1 的列向量，空间网格剖分 % T 长度为 NT +1 的行向量，时间网格剖分
% U (NS +1) *( NT +1) 矩阵，U(:,i)表示第i个时间层网格部分上的数值解 [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
display('时间空间离散不满足向前差分的稳定条件'); end
U=zeros(N,M);
U(:,1)=pde.u_initial(X); U(1,:)=pde.u_left(T); U(end,:)=pde.u_right(T); %% 向前差分
d=1-2*ones(1,N-2)*r; c=ones(1,N-3)*r;
A=diag(c,-1)+diag(d)+diag(c,1); for i=2:M
RHS=tau*pde.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 emax=heat_equation_fd1d_error(X,T,U,u_exact)
ue=u_exact(X,T); %真解在网格节点处的值 emax=max(max(ue-U)); end
function showsolution (X,T,U) %% SHOWSOLUTION 以二元函数方式显示数值解 % 输入参数
% X 长度为 NS +1 的列向量，空间网格剖分N % T 长度为 NT +1 的行向量，时间网格剖分M
% U N*M 矩阵，U(:,i)表示第i个时间层网格部分上的数值解 [x,t] = meshgrid (X,T); mesh (x,t,U'); xlabel ('X'); ylabel ('T'); zlabel ('U(X,T)'); End
function showvarysolution (X,T,U) %% SHOWVARYSOLUTION 显示数值解随着时间的变化 % % 输入参数
% X 长度为 NS +1 的列向量，空间网格剖分N % T 长度为 NT +1 的行向量，时间网格剖分M
% U N*M 矩阵，U(:,i)表示第i个时间层网格部分上的数值解 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 u_exact=u_exact(X,T) u_exact=sin(2*pi*X)*exp(10*T); End
function pde=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);%T为行向量 tau=0.1/(NT-1); end
function [X,h]=space_grid(NS) X=linspace(0,1,NS)';%X为列向量 h=1/(NS-1); end
function u=u_initial(X) u=sin(2*pi*X); end
function u=u_left(T)%左边值条件，x0位置不动，时间在变化 u=zeros(size(T)); end
function u=u_right(T)%右边值条件，xN位置不动，时间在变化 u=zeros(size(T)); end
function f=f(x,t)
f=10*exp(10*t)*sin(2*pi*x)+4*a*pi*pi*exp(10*t)*sin(2*pi*x); end
function a=a() a=1; end end
pde=model_data(); %模型数据结构体 main_test %向前差分
N=[11,21,31]; M=[101,201,301]; for i=1:3
[X,T,U]=heat_equation_fd1d(N(i),M(i),pde); showvarysolution(X,T,U); %以随时间变化方式显示数值解 showsolution(X,T,U);%以二元函数方式显示数值解
[emax]=heat_equation_fd1d_error(X,T,U,@u_exact); format short disp('emax'); disp([emax]); end
实验结果：在命令窗口输入>> main_test 回车可得运算结果和图像 >> main_test emax