一维热传导方程一． 问题介绍考虑一维热传导方程： （1）,0),(22T t x f xu atu ≤<+∂∂=∂∂其中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 都是正整数。

用两族平行直线：),,1,0(N j jh x x j ===),,1,0(M k k t t k ===τ将矩形域}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 kj k j k j kjk jf hu u u au u ++-=--++τ)(j j x f f =，)(0j j j x u ϕϕ==， 00==kN k u u ，其中j = 1,2,…，N-1，k = 1,2,…,M-1。

以2/h a r τ=表示网比。

则方程（5）可以改写为：j k j kj kj k jf ruu r ruu τ++-+=-++111)21(易知向前差分格式是显格式。

2. 向后差分格式 （6） ,11111)21(j kj k j k jk j f u ruu u ruτ+=-++-+-+++)(0j j j x u ϕϕ==， 00==kN ku u ，其中j = 1,2,…，N-1，k = 1,2,…,M-1，易知向前差分格式是显格式。

3. 六点对称格式（Grank-Nicolson 格式）将向前差分格式和向后差分格式作算术平均，即得到六点对称格式：（7） 111112)1(2+-+++-++-k j k jk j ur ur ur=j k j kj kj f ur u r u rτ++-+-+112)1(2利用0j u 和边值便可逐层求到kj u 。

六点对称格式是隐格式，由第k 层计算第k+1层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。

将其截断误差)(x R kj 于),(21+k j tx （21+k t=τ)21(+k ）展开，则得)(x R k j =)(22h O +τ。

4. Richardson 格式 （8） τ211-+-k j k ju u j kj k j k j f hu u u a++-=-+2112，或（9）j k jk j k j k j k jf u u u u r u τ2)2(21111+++-=--++。

这是三层显示差分格式。

截断误差阶为)(22h O +τ。

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

四．格式稳定性通过误差估计方程1111)2(2--++++-=k jkj kj kj k je e e e r e（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] 0.333333 0.055556 1.475341 1.473867 0.001474 0.666667 0.055556 2.059004 2.056947 0.002057 0.333333 0.111111 1.559623 1.557037 0.002586 0.666667 0.111111 2.176630 2.173719 0.002911 0.333333 0.166667 1.648721 1.645619 0.003102 0.666667 0.166667 2.300976 2.297385 0.003591 0.333333 0.222222 1.742909 1.739373 0.003536 0.666667 0.222222 2.432425 2.428445 0.003981 0.333333 0.277778 1.842477 1.838647 0.003831 0.666667 0.277778 2.571384 2.567048 0.004337 0.333333 0.333333 1.947734 1.943620 0.004114 0.666667 0.333333 2.718282 2.713651 0.004630 0.333333 0.388889 2.059004 2.054632 0.004372 0.666667 0.388889 2.873571 2.868644 0.004927 0.333333 0.444444 2.176630 2.171992 0.004638 0.666667 0.444444 3.037732 3.032512 0.0052200.333333 0.500000 2.300976 2.296068 0.004908 0.666667 0.500000 3.211271 3.205744 0.005526 0.333333 0.555556 2.432425 2.427233 0.005193 0.666667 0.555556 3.394723 3.388878 0.005845 0.333333 0.611111 2.571384 2.565894 0.005491 0.666667 0.611111 3.588656 3.582475 0.006181 0.333333 0.666667 2.718282 2.712476 0.005805 0.666667 0.666667 3.793668 3.787133 0.006535 0.333333 0.722222 2.873571 2.867434 0.006137 0.666667 0.722222 4.010392 4.003483 0.006908 0.333333 0.777778 3.037732 3.031243 0.006488 0.666667 0.777778 4.239496 4.232193 0.007303 0.333333 0.833333 3.211271 3.204411 0.006859 0.666667 0.833333 4.481689 4.473969 0.007721 0.333333 0.888889 3.394723 3.387472 0.007251 0.666667 0.888889 4.737718 4.729556 0.008162 0.333333 0.944444 3.588656 3.580991 0.007665 0.666667 0.944444 5.008373 4.999745 0.008628 当n等于2和m等于17时最大误差为0.008628其中r = 1/2，格式是稳定的。

2. 用向后差分格式验证得数值结果如下：请输入n的值(输入0结束程序):6请输入m的值(输入0结束程序):6xj 真实值x[i] 近似值u[i] 误差err[i]第1层结果时间节点Tk=0.1428570.142857 1.330712 1.335002 0.0042890.285714 1.535063 1.542358 0.0072950.428571 1.770795 1.779949 0.0091540.571429 2.042727 2.052524 0.0097970.714286 2.356418 2.365346 0.0089270.857143 2.718282 2.724256 0.005974第1层的最大误差是0.009797第2层结果时间节点Tk=0.2857140.142857 1.535063 1.541797 0.0067340.285714 1.770795 1.782424 0.0116290.428571 2.042727 2.057345 0.0146180.571429 2.356418 2.371895 0.0154770.714286 2.718282 2.732070 0.0137880.857143 3.135715 3.144633 0.008919第2层的最大误差是0.015477第3层结果时间节点Tk=0.4285710.285714 2.042727 2.057518 0.014791 0.428571 2.356418 2.375010 0.018592 0.571429 2.718282 2.737884 0.019602 0.714286 3.135715 3.153041 0.017326 0.857143 3.617251 3.628337 0.011086 第3层的最大误差是0.019602第4层结果时间节点Tk=0.571429 0.142857 2.042727 2.052888 0.010161 0.285714 2.356418 2.374062 0.017644 0.428571 2.718282 2.740457 0.022175 0.571429 3.135715 3.159058 0.023343 0.714286 3.617251 3.637827 0.020576 0.857143 4.172734 4.185851 0.013117 第4层的最大误差是0.023343第5层结果时间节点Tk=0.714286 0.142857 2.356418 2.368276 0.011857 0.285714 2.718282 2.738880 0.020598 0.428571 3.135715 3.161600 0.025886 0.571429 3.617251 3.644485 0.027234 0.714286 4.172734 4.196716 0.023982 0.857143 4.813520 4.828788 0.015268 第5层的最大误差是0.027234第6层结果时间节点Tk=0.857143 0.142857 2.718282 2.732017 0.013735 0.285714 3.135715 3.159578 0.023864 0.428571 3.617251 3.647240 0.029989 0.571429 4.172734 4.204278 0.031544 0.714286 4.813520 4.841287 0.027767 0.857143 5.552708 5.570378 0.017670 第6层的最大误差是0.031544当n等于6时最大误差为0.0315443. 用六点对称格式验证数值结果如下：请输入n的值(输入0结束程序):6请输入m的值(输入0结束程序):6xj 真实值x[i] 近似值u[i] 误差err[i] 第1层结果时间节点Tk=0.142857 0.142857 1.330712 1.330988 0.000276 0.285714 1.535063 1.535522 0.000459 0.428571 1.770795 1.771368 0.000573 0.571429 2.042727 2.043350 0.0006230.857143 2.718282 2.718692 0.000410 第1层的最大误差是0.000623第2层结果时间节点Tk=0.285714 0.142857 1.535063 1.535424 0.000361 0.285714 1.770795 1.771435 0.000640 0.428571 2.042727 2.043538 0.000810 0.571429 2.356418 2.357268 0.000849 0.714286 2.718282 2.719016 0.000734 0.857143 3.135715 3.136162 0.000447 第2层的最大误差是0.000849第3层结果时间节点Tk=0.428571 0.142857 1.770795 1.771233 0.000438 0.285714 2.042727 2.043476 0.000749 0.428571 2.356418 2.357355 0.000937 0.571429 2.718282 2.719268 0.000986 0.714286 3.135715 3.136592 0.000877 0.857143 3.617251 3.617825 0.000574 第3层的最大误差是0.000986第4层结果时间节点Tk=0.571429 0.142857 2.042727 2.043222 0.000495 0.285714 2.356418 2.357287 0.000869 0.428571 2.718282 2.719377 0.001095 0.571429 3.135715 3.136867 0.001153 0.714286 3.617251 3.618261 0.001010 0.857143 4.172734 4.173365 0.000631 第4层的最大误差是0.001153第5层结果时间节点Tk=0.714286 0.142857 2.356418 2.356999 0.000581 0.285714 2.718282 2.719284 0.001003 0.428571 3.135715 3.136973 0.001258 0.571429 3.617251 3.618573 0.001322 0.714286 4.172734 4.173900 0.001166 0.857143 4.813520 4.814271 0.000751 第5层的最大误差是0.001322第6层结果时间节点Tk=0.857143 0.142857 2.718282 2.718945 0.000663 0.285714 3.135715 3.136871 0.001157 0.428571 3.617251 3.618705 0.001455 0.571429 4.172734 4.174265 0.001531 0.714286 4.813520 4.814867 0.001347 0.857143 5.552708 5.553559 0.000851 第6层的最大误差是0.001531当n等于6时最大误差为0.0015314．用Richardson格式验证数值结果如下：请输入n的值(输入0结束程序):5请输入m的值(输入0结束程序):5xj tk 真实值x[i][k] 近似值u[i][k] 误差err[i][k]0.166667 0.166667 1.395612 1.396080 0.0004680.333333 0.166667 1.648721 1.649481 0.0007600.500000 0.166667 1.947734 1.948649 0.0009150.666667 0.166667 2.300976 2.301888 0.0009120.833333 0.166667 2.718282 2.718949 0.0006670.166667 0.333333 1.648721 1.645540 0.0031820.333333 0.333333 1.947734 1.944806 0.0029280.500000 0.333333 2.300976 2.297583 0.0033930.666667 0.333333 2.718282 2.713600 0.0046820.833333 0.333333 3.211271 3.204099 0.0071710.166667 0.500000 1.947734 1.988145 0.0404110.333333 0.500000 2.300976 2.291621 0.0093550.500000 0.500000 2.718282 2.707514 0.0107680.666667 0.500000 3.211271 3.195685 0.0155860.833333 0.500000 3.793668 3.907779 0.1141110.166667 0.666667 2.300976 1.214156 1.0868200.333333 0.666667 2.718282 3.293822 0.5755410.500000 0.666667 3.211271 3.164907 0.0463640.666667 0.666667 3.793668 5.400692 1.6070240.833333 0.666667 4.481689 1.545878 2.9358110.166667 0.833333 2.718282 35.747074 33.0287920.333333 0.833333 3.211271 -24.211361 27.4226310.500000 0.833333 3.793668 31.083927 27.2902590.666667 0.833333 4.481689 -69.891509 74.3731980.833333 0.833333 5.294490 95.148891 89.854401附录1#include <stdio.h>#include <stdlib.h>#include <math.h>#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; //将空间区间【0，1】分为n等份；时间区间【0，1】分为m等份void catchup(){int i;beta[1]=c[1]/b[1];for(i=2;i<n;i++)beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);y[1]=f[1]/b[1];for(i=2;i<=n;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);u[n][1]=y[n];for(i=n-1;i>0;i--)u[i][1]=y[i]-beta[i]*u[i+1][1];}int main() //一维热传导方程的向前差分格式{int k,i;double h,t,r;double pi=3.1415627;printf("请输入n的值(输入0结束程序):\n");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):\n");while(scanf("%d",&m)&&m&&n) //u(x,t)=exp(x+t),u(x,0)=exp(x),f(x)=0,x属于[0,1],t属于[0,1],a=1.{h=1.0/(n+1);t=1.0/(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");for(k=1;k<=m;k++){for(int j=1;j<=n;j++){u[j][k]=r*u[j+1][k-1]+(1-2*r)*u[j][k-1]+r*u[j-1][k-1];}}double T=0;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]);if(err[i][k]>T) T=err[i][k];}}printf("当n等于%d和m等于%d时最大误差为%lf\n",n,m,T);printf("请输入n的值(输入0结束程序):");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):");}system("PASUE");return 0;}附录2#include <stdio.h>#include <stdlib.h>#include <math.h>#define Max_N 1000double u[Max_N],b[Max_N],a[Max_N],c[Max_N],f[Max_N],err[Max_N],x[Max_N],y[Max_N],beta[Max_N],Err[Max_N];int n,m; //将空间区间【0，1】分为n等份；时间区间【0，1】分为m等份void catchup(){int i;beta[1]=c[1]/b[1];for(i=2;i<n;i++)beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);y[1]=f[1]/b[1];for(i=2;i<=n;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);u[n]=y[n];for(i=n-1;i>0;i--)u[i]=y[i]-beta[i]*u[i+1];}int main() //一维热传导方程的六点对称格式{int k,i;double h,t,r;//int j=0;double pi=3.1415627;printf("请输入n的值(输入0结束程序):\n");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):\n");while(scanf("%d",&m)&&m&&n) //u(x,t)=exp(x+t),u(x,0)=exp(x),f(x)=0,x属于[0,1],t属于[0,1],a=1.{h=1.0/(n+1);t=1.0/(m+1);r=t/(h*h);printf("xj 真实值x[i] 近似值u[i] 误差err[i]\n");double M=0;//double temp=0;for(k=1;k<=m;k++){b[1]=1+2*r;c[1]=-r;a[n]=-r;b[n]=1+2*r;if(k==1){f[1]=exp(h)+r*exp(t);f[n]=exp(n*h)+r*exp((n+1)*h+t);}else{f[1]=u[1]+r*exp(k*t);f[n]=u[n]+r*exp((n+1)*h+k*t);}for(i=2;i<n;i++){b[i]=1+2*r;a[i]=-r;c[i]=-r;if(k==1) f[i]=exp(i*h);else f[i]=u[i];}catchup();printf("第%d层结果时间节点Tk=%lf\n",k,k*t);double T=0;//double tem=0;for(i=1;i<=n;i++){x[i]=exp(i*h+k*t);err[i]=x[i]>u[i]?x[i]-u[i]:u[i]-x[i];printf("%lf %lf %lf %lf\n",i*h,x[i],u[i],err[i]);if(err[i]>T) T=err[i];//tem=tem+err[i]*err[i]*h;}printf("第%d层的最大误差是%lf\n",k,T);if(T>M) M=T;//temp=temp+tem;}printf("\n");//Err[j]=temp;printf("当n等于%d时最大误差为%lf\n",n,M);//printf("当n等于%d时其二范数为%lf\n",n,Err[j]);//if(j!=0) printf("log(Err[j-1]/Err[j])=%lf\n\n",log(Err[j-1]/Err[j]));//j++;printf("请输入n的值(输入0结束程序):");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):");}system("PASUE");return 0;}附录3#include <stdio.h>#include <stdlib.h>#include <math.h>#define Max_N 1000double u[Max_N],b[Max_N],a[Max_N],c[Max_N],f[Max_N],err[Max_N], x[Max_N],y[Max_N],beta[Max_N],Err[Max_N];int n,m; //将空间区间【0，1】分为n等份；时间区间【0，1】分为m等份void catchup(){int i;beta[1]=c[1]/b[1];for(i=2;i<n;i++)beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);y[1]=f[1]/b[1];for(i=2;i<=n;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);u[n]=y[n];for(i=n-1;i>0;i--)u[i]=y[i]-beta[i]*u[i+1];}int main() //一维热传导方程的六点对称格式{int k,i;double h,t,r;//int j=0;double pi=3.1415627;printf("请输入n的值(输入0结束程序):\n");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):\n");while(scanf("%d",&m)&&m&&n) //u(x,t)=exp(x+t),u(x,0)=exp(x),f(x)=0,x属于[0,1],t属于[0,1],a=1.{h=1.0/(n+1);t=1.0/(m+1);r=t/(h*h);printf("xj 真实值x[i] 近似值u[i] 误差err[i]\n");double M=0;//double temp=0;for(k=1;k<=m;k++){b[1]=1+r;c[1]=-r/2;a[n]=-r/2;b[n]=1+r;if(k==1){f[1]=r/2*exp(2*h)+(1-r)*exp(h)+r/2+r/2*exp(t);f[n]=r/2*exp((n+1)*h)+(1-r)*exp(n*h)+r/2*exp((n-1)*h)+r/2*exp((n+1)*h+t);}else{f[1]=r/2*u[2]+(1-r)*u[1]+r/2*exp((k-1)*t)+r/2*exp(k*t);f[n]=r/2*u[n-1]+(1-r)*u[n]+r/2*exp((n+1)*h+(k-1)*t)+r/2*exp((n+1)*h+k*t);}for(i=2;i<n;i++){b[i]=1+r;a[i]=-r/2;c[i]=-r/2;if(k==1) f[i]=r/2*exp((i+1)*h)+(1-r)*exp(i*h)+r/2*exp((i-1)*h);else f[i]=r/2*u[i+1]+(1-r)*u[i]+r/2*u[i-1];}catchup();printf("第%d层结果时间节点Tk=%lf\n",k,k*t);double T=0;double tem=0;for(i=1;i<=n;i++){x[i]=exp(i*h+k*t);err[i]=x[i]>u[i]?x[i]-u[i]:u[i]-x[i];printf("%lf %lf %lf %lf\n",i*h,x[i],u[i],err[i]);if(err[i]>T) T=err[i];//tem=tem+err[i]*err[i]*h;}printf("第%d层的最大误差是%lf\n",k,T);if(T>M) M=T;//temp=temp+tem;}printf("\n");//Err[j]=temp;printf("当n等于%d时最大误差为%lf\n",n,M);//printf("当n等于%d时其二范数为%lf\n",n,Err[j]);//if(j!=0) printf("log(Err[j-1]/Err[j])=%lf\n\n",log(Err[j-1]/Err[j]));//j++;printf("请输入n的值(输入0结束程序):");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):");}system("PASUE");return 0;}附录4#include <stdio.h>#include <stdlib.h>#include <math.h>#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; //将空间区间【0，1】分为n等份；时间区间【0，1】分为m等份void catchup(){int i;beta[1]=c[1]/b[1];for(i=2;i<n;i++)beta[i]=c[i]/(b[i]-a[i]*beta[i-1]);y[1]=f[1]/b[1];for(i=2;i<=n;i++)y[i]=(f[i]-a[i]*y[i-1])/(b[i]-a[i]*beta[i-1]);u[n][1]=y[n];for(i=n-1;i>0;i--)u[i][1]=y[i]-beta[i]*u[i+1][1];}int main() //一维热传导方程的Richardson格式{int k,i;double h,t,r;double pi=3.1415627;printf("请输入n的值(输入0结束程序):\n");if(scanf("%d",&n)) printf("请输入m的值(输入0结束程序):\n");while(scanf("%d",&m)&&m&&n) //u(x,t)=exp(x+t),u(x,0)=exp(x),f(x)=0,x属于[0,1],t属于[0,1],a=1.{h=1.0/(n+1);t=1.0/(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;}。