一维热传导方程一． 问题介绍考虑一维热传导方程：（1） ,0),(22T t x f xu 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 的网点集合；h Γ=h G --h G 是网格界点集合。

三． 离散格式第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。

第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。

1. 向前差分格式（5） ,22111j k j k j kj kjk j f h u u u a u u ++-=--++τ)(j j x f f =， )(0j 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 τ+=-++-+-+++ )(0j 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 kj 于),(21+k j tx （21+k t =τ)21(+k ）展开，则得 )(x R kj =)(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 +τ。

为了使计算能够逐层进行，除初值0j u 外，还要用到1j u ，这可以用前述二层差分格式计算（为保证精度，可将[0,τ]分成若干等份）。

四． 格式稳定性通过误差估计方程（1） 可知对任意的r ，Richardson 格式都不稳定，所以Richardson 格式绝对不稳定。

（2） 当210≤<r 时，向前差分格式趋于稳定；当21>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结束程序):17xj tk 真实值x[i][k]近似值u[i][k]误差err[i][k]当n等于2和m等于17时最大误差为其中r = 1/2，格式是稳定的。

2. 用向后差分格式验证得数值结果如下：请输入n的值(输入0结束程序):6请输入m的值(输入0结束程序):6xj 真实值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结束程序):6xj 真实值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结束程序):5xj tk 真实值x[i][k] 近似值u[i][k] 误差err[i][k]附录1#include <>#include <>#include <>#define Max_N 1000double 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<n;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;}。