一维热传导方程
一维热传导方程Last revision on 21 December 2020
一维热传导方程
一． 问题介绍
考虑一维热传导方程：
（1） ,0),(22T t x f x
u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类：
第一类、初值问题（也称Cauthy 问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（∞<<∞-x ）和初始条件：
（2）
),()0,(x x u ?= ∞<<∞-x 第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数),(t x u ，满足方程（1）（l x <<0）和初始条件：
（3）
),()0,(x x u ?= l x <<0 及边值条件
(4)
.0),(),0(==t l u t u T t ≤≤0 假定)(x ?在相应区域光滑，并且在l x ,0=满足相容条件，使上述问题有唯一充分光滑的解。
二． 区域剖分
考虑边值问题（1），（4）的差分逼近。去空间步长N l h /=和时间步长M T /=τ，其中N,M 都是正整数。用两族平行直线： 将矩形域}0;0{T t l x G ≤≤≤≤=分割成矩形网格，网格节点为),(k j t x 。以h G 表示网格内点集合，即位于开矩形G 的网点集合；h G 表示所有位于闭矩形G 的网点集合；Γ=G --G 是网格界点集合。
三． 离散格式
第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。
第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。
1. 向前差分格式
（5）
,22111j k j k j k j k j k j f h u u u a u u ++-=--++τ
)(j j x f f =， )(0
j j j x u ??==， 00==k N k u u ，
其中j = 1,2,…，N-1，k = 1,2,…,M-1。以2/h a r τ=表示网比。则方程（5）可以改写为：
易知向前差分格式是显格式。
2. 向后差分格式
（6）
,11111)21(j k j k j k j k j f u ru u u ru τ+=-++-+-+++ )(0
j j j x u ??==， 00==k N k u u ，
其中j = 1,2,…，N-1，k = 1,2,…,M-1，易知向前差分格式是显格式。
3. 六点对称格式（Grank-Nicolson 格式）
将向前差分格式和向后差分格式作算术平均，即得到六点对称格式：
（7）
111112)1(2+-+++-++-k j k j k j u r u r u r =j k j k j k
j f u r u r u r τ++-+-+112
)1(2 利用0j
u 和边值便可逐层求到k j u 。六点对称格式是隐格式，由第k 层计算第k+1层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。
将其截断误差)(x R k
j 于),(21+k j t
x （21+k t =τ)21(+k ）展开，则得 )(x R k
j =)(22h O +τ 。
4. Richardson 格式
（8）
τ211-+-k j k j u u j k j k j k j f h u u u a ++-=-+2112，
或
（9） j k j k j k j k j k j f u u u u r u τ2)2(21111+++-=--++。 这是三层显示差分格式。截断误差阶为)(22h O +τ。为了使计算能够逐层进
行，除初值0
j u 外，还要用到1j u ，这可以用前述二层差分格式计算（为保证精
度，可将[0,τ]分成若干等份）。
四． 格式稳定性
通过误差估计方程
（1） 可知对任意的r ，Richardson 格式都不稳定，所以Richardson 格式绝对不稳定。
（2） 当210≤ 1>r 时，向前差分格式的误差无限增长。因此向前差分格式是条件稳定。 （3） 向后差分格式和六点对称格式都绝对稳定，且各自的截断误差阶分别为)(2h O +τ和)(22h O +τ。 五． 数值例子 例1 令f ( x ) = 0和a = 1，可求得u （x,t ）一个解析解为u ( x , t ) =exp ( x + t )。 1．用向前差分格式验证得数值结果如下： 请输入n 的值(输入0结束程序): 2 请输入m的值(输入0结束程序): 17 xj tk 真实值x[i][k] 近似值u[i][k] 误差err[i][k] 当n等于2和m等于17时最大误差为 其中r = 1/2，格式是稳定的。 2. 用向后差分格式验证得数值结果如下：请输入n的值(输入0结束程序): 6 请输入m的值(输入0结束程序): 6 xj 真实值x[i] 近似值u[i] 误差err[i] 第1层结果时间节点Tk= 第1层的最大误差是 第2层结果时间节点Tk= 第2层的最大误差是 第3层结果时间节点Tk= 第3层的最大误差是 第4层结果时间节点Tk= 第4层的最大误差是 第5层结果时间节点Tk= 第5层的最大误差是 第6层结果时间节点Tk= 第6层的最大误差是 当n等于6时最大误差为 3. 用六点对称格式验证数值结果如下： 请输入n的值(输入0结束程序): 6 请输入m的值(输入0结束程序): 6 xj 真实值x[i] 近似值u[i] 误差err[i] 第1层结果时间节点Tk= 第1层的最大误差是 第2层结果时间节点Tk= 第2层的最大误差是 第3层结果时间节点Tk= 第3层的最大误差是 第4层结果时间节点Tk= 第4层的最大误差是 第5层结果时间节点Tk= 第5层的最大误差是 第6层结果时间节点Tk= 第6层的最大误差是 当n等于6时最大误差为 4．用Richardson格式验证数值结果如下： 请输入n的值(输入0结束程序): 5 请输入m的值(输入0结束程序): 5 xj tk 真实值x[i][k] 近似值u[i][k] 误差err[i][k] 附录1 #include <> #include <> #include <> #define Max_N 1000 double u[Max_N][Max_N],b[Max_N],a[Max_N],c[Max_N],f[Max_N], err[Max_N][Max_N],x[Max_N][Max_N],y[Max_N],beta[Max_N],Err[Max_ N]; int n,m; { h=(n+1); t=(m+1); r=t/(h*h); for(i=0;i<=n+1;i++) { h=(n+1); t=(m+1); r=t/(h*h); printf("xj 真实值x[i] 近似值u[i] 误差err[i]\n"); double M=0; { h=(n+1); t=(m+1); r=t/(h*h); printf("xj 真实值x[i] 近似值u[i] 误差err[i]\n"); double M=0; { h=(n+1); t=(m+1); r=t/(h*h); for(i=0;i<=n+1;i++)//初值条件 { u[i][0]=exp(i*h); } for(k=0;k<=m+1;k++)//边值条件 { u[0][k]=exp(k*t); u[n+1][k]=exp((n+1)*h+k*t); } printf("xj tk 真实值x[i][k] 近似值u[i][k] 误差 err[i][k]\n"); b[1]=1+r; c[1]=-r/2; a[n]=-r/2; b[n]=1+r; f[1]=r/2*u[2][0]+(1-r)*u[1][0]+r/2+r/2*u[0][1]; f[n]=r/2*u[n+1][0]+(1-r)*u[n][0]+r/2*u[n-1][0]+r/2*u[n+1][1]; for(i=2;i { b[i]=1+r; a[i]=-r/2; c[i]=-r/2; f[i]=r/2*u[i+1][0]+(1-r)*u[i][0]+r/2*u[i-1][0]; } catchup(); for(k=2;k<=m;k++) { for(int j=1;j<=n;j++) { u[j][k]=2*r*(u[j+1][k-1]-2*u[j][k-1]+u[j-1][k-1])+u[j][k-2]; } } for(k=1;k<=m;k++) { for(i=1;i<=n;i++) { x[i][k]=exp(i*h+k*t); err[i][k]=x[i][k]>u[i][k]x[i][k]-u[i][k]:u[i][k]-x[i][k]; printf("%lf %lf %lf%lf %lf\n",i*h,k*t,x[i][k],u[i][k],err[i][k]); } } printf("请输入n的值(输入0结束程序):"); if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):"); } system("PASUE"); return 0; }