一维热传导方程(Richardson格式)
中南林业科技大学
偏微分方程数值解法学生姓名：周晓虹
学号：20083710
学院：理学院
专业年级：08信计1班
设计题目：一维热传导方程的Richardson格式
2011年06月
一． 问题介绍
考虑一维热传导方程： （1）
,0),(22
T 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 都是正整数。用两族平行直线：
),,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 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。
Richardson 格式 （8） τ
211-+-k j
k j
u u
j k j k j
k j f h
u
u u
a
++-=-+2
1
1
2，
或 （9）
j k j
k
j k j
k j k j
f u u
u
u
r u
τ2)2(211
1
1+++-=--++。
这是三层显示差分格式。
截断误差阶为：
)(2
2
h O +τ
。
为了使计算能够逐层进行，除初值0
j u 外，还要用到1
j u ，这可以用前述二层差分格式计算（为保证精度，可将[0,τ]分成若干等份）。
四． 格式稳定性
通过误差估计方程
1
111
)2(2--++++-=k j
k j k j k j k j
e e e e r e
（1） 可知对任意的r ，Richardson 格式都不稳定，所以Richardson 格式绝对不稳定。
五． 数值例子
例1 令f ( x ) = 0和a = 1，可求得u （x,t ）一个解析解为u ( x , t ) =exp ( x + t )。 用Richardson 格式验证数值结果如下： 请输入n 的值(输入0结束程序): 5
请输入m 的值(输入0结束程序): 5
xj tk 真实值x[i][k] 近似值u[i][k] 误差err[i][k] 0.166667 0.166667 1.395612 1.396080 0.000468 0.333333 0.166667 1.648721 1.649481 0.000760 0.500000 0.166667 1.947734 1.948649 0.000915 0.666667 0.166667 2.300976 2.301888 0.000912 0.833333 0.166667 2.718282 2.718949 0.000667 0.166667 0.333333 1.648721 1.645540 0.003182
0.333333 0.333333 1.947734 1.944806 0.002928
0.500000 0.333333 2.300976 2.297583 0.003393
0.666667 0.333333 2.718282 2.713600 0.004682
0.833333 0.333333 3.211271 3.204099 0.007171
0.166667 0.500000 1.947734 1.988145 0.040411
0.333333 0.500000 2.300976 2.291621 0.009355
0.500000 0.500000 2.718282 2.707514 0.010768
0.666667 0.500000 3.211271 3.195685 0.015586
0.833333 0.500000 3.793668 3.907779 0.114111
0.166667 0.666667 2.300976 1.214156 1.086820
0.333333 0.666667 2.718282 3.293822 0.575541
0.500000 0.666667 3.211271 3.164907 0.046364
0.666667 0.666667 3.793668 5.400692 1.607024
0.833333 0.666667 4.481689 1.545878 2.935811
0.166667 0.833333 2.718282 35.747074 33.028792
0.333333 0.833333 3.211271 -24.211361 27.422631
0.500000 0.833333 3.793668 31.083927 27.290259
0.666667 0.833333 4.481689 -69.891509 74.373198
0.833333 0.833333 5.294490 95.148891 89.854401
六．参考文献
[1]陆金甫，关治编.偏微分方程数值解法[M].北京：清华大学出版社，2003：35-137.
[2]南京大学数学系.计算数学专业编，偏微分方程数值解法[M].北京：科学出版社，1979：10-11.
附录
程序源代码：
#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; //将空间区间【0，1】分为n等份；时间区间【0，1】分为m等份
void catchup()
{
int i;
beta[1]=c[1]/b[1];
for(i=2;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 { 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; }