一维热传导方程

一维热传导方程

一． 问题介绍

考虑一维热传导方程： （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 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。

1. 向前差分格式 （5） ,22

1

11

j 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）可以改写为：

j k j k

j k

j k j

f ru

u r ru

u τ++-+=-++1

1

1

)21(

易知向前差分格式是显格式。

2. 向后差分格式 （6） ,11

1

11

)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） 11

111

2

)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

τ++

-+-+1

12

)1(2

利用0

j u 和边值便可逐层求到k

j u 。六点对称格式是隐格式，由第k 层计算第k+1层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。 将其截断误差)(x R k

j 于),(2

1+

k j t

x （2

1

+

k t

=τ)2

1(+k ）展开，则得

)(x R k j =)(2

2

h O +τ

。

4. Richardson 格式 （8） τ

21

1

-+-k j k j

u u j k

j k j k j f h

u u u a

++-=-+2

1

12，

或

（9）

j k j

k j k j k j k j

f u u u u r u τ2)2(21

111

+++-=--++。

这是三层显示差分格式。截断误差阶为)(2

2

h O +τ

。为了使计算能够逐层进行，除初值0j

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 格式绝对不稳定。 （2） 当2

10≤

1>

r 时，向前差分格式的误差无限增

长。因此向前差分格式是条件稳定。

（3） 向后差分格式和六点对称格式都绝对稳定，且各自的截断误差阶分别为)(2h O +τ和

)(2

2

h 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] 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.005220

0.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结束程序):

6

xj 真实值x[i] 近似值u[i] 误差err[i]

第1层结果时间节点Tk=0.142857

0.142857 1.330712 1.335002 0.004289

0.285714 1.535063 1.542358 0.007295

0.428571 1.770795 1.779949 0.009154

0.571429 2.042727 2.052524 0.009797

0.714286 2.356418 2.365346 0.008927

0.857143 2.718282 2.724256 0.005974

第1层的最大误差是0.009797

第2层结果时间节点Tk=0.285714

0.142857 1.535063 1.541797 0.006734

0.285714 1.770795 1.782424 0.011629

0.428571 2.042727 2.057345 0.014618

0.571429 2.356418 2.371895 0.015477

0.714286 2.718282 2.732070 0.013788

0.857143 3.135715 3.144633 0.008919

第2层的最大误差是0.015477

第3层结果时间节点Tk=0.428571

0.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.031544

3. 用六点对称格式验证数值结果如下：

请输入n的值(输入0结束程序):

6

请输入m的值(输入0结束程序):

6

xj 真实值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.000623

0.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.001531

4．用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

#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() //一维热传导方程的向前差分格式

{

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

#include

#include

#define Max_N 1000

double 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

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

{

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

#include

#include

#define Max_N 1000

double 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

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

{

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

#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;

}

一维热传导方程的差分格式

《微分方程数值解》 课程论文 学生姓名1：许慧卿学号：20144329 学生姓名2：向裕学号：20144327学生姓名3：邱文林学号：20144349学生姓名4：高俊学号：20144305学生姓名5：赵禹恒学号：20144359学生姓名6：刘志刚学号： 20144346 学院：理学院 专业：14级信息与计算科学 指导教师：陈红斌 2017年6 月25日

《偏微分方程数值解》课程论文 《一维热传导方程的差分格式》论文 一、《微分方程数值解》课程论文的格式 1）引言：介绍研究问题的意义和现状 2）格式：给出数值格式 3）截断误差：给出数值格式的截断误差 4）数值例子：按所给数值格式给出数值例子 5）参考文献：论文所涉及的文献和教材 二、《微分方程数值解》课程论文的评分标准 1）文献综述：10分； 2）课题研究方案可行性：10分； 3）数值格式：20分； 4）数值格式的算法、流程图：10分； 5）数值格式的程序：10分； 6）论文撰写的条理性和完整性：10分； 7）论文工作量的大小及课题的难度：10分； 8）课程设计态度：10分； 9）独立性和创新性：10分。 评阅人： - 2 -

一维热传导方程的差分格式 1 引言 考虑如下一维非齐次热传导方程Dirichlet 初边值问题 22(,),u u a f x t t x ??=+?? ,c x d << 0,t T <≤ (1.1) (,0)(),u x x ?= ,c x d ≤≤ (1.2) (,)(),u c t t α= (,)(),u d t t β= 0t T <≤ (1.3) 的有限差分方法, 其中a 为正常数,(,),(),(), ()f x t x t t ?αβ为已知常数, ()(0),c ?α= ()(0).d ?β= 称(1.2)为初值条件, (1.3)为边值条件. 本文将给出(1.1) (1.3)的向前Euler 格式, 向后Euler 格式和Crank Nicolson -格式, 并给出其截断误差和数值例子. 经对比发现, Crank Nicolson -格式误差最小, 向前 Euler 格式次之, 向后Euler 格式误差最大. 2 差分格式的建立 2.1 向前Euler 格式 将区间[,]c d 作M 等分, 将[]0,T 作N 等分, 并记 ()/h d c M =-, /T N τ=, j x c jh =+,0j M ≤≤, k t k τ=,0k N ≤≤. 分别称h 和τ为空间步长和时间步长.用 两组平行直线 j x x =, 0j M ≤≤, k t t =, 0k N ≤≤ 将Ω分割成矩形网格.记{} |0h j x j M Ω=≤≤, {}|0k t k N τΩ=≤≤, h h ττΩ=Ω?Ω. 称() ,j k x t 为结点[1] . 定义h τΩ上的网格函数 {}|0,0k j U j M k N Ω=≤≤≤≤, 其中() ,k j j k U u x t =. 在结点() ,j k x t 处考虑方程(1.1),有

一维热传导方程

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （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 的网点集合；h Γ=h G --h G 是网格界点集合。 三． 离散格式 第k+1层值通过第k 层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。 第k+1层值不能通过第k 层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。 1. 向前差分格式 （5） ,221 11j k j k j k j k j k j f h u u u a u u ++-=--++τ

一维热传导方程

一维热传导方程 一． 问题介绍 考虑一维热传导方程： （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 层值明显表示出来，无需求解线性代数方程组，这样的格式称为

偏微分一维热传导问题

偏微分大作业 一维热传导方程问题——运用隐式格式求解数值解

目录 问题描述 (3) 1 解析解——分离变量法 (3) 2 数值解——隐式格式 (5) 3 证明隐式格式的相容性与稳定性 (5) 4 数值解——分析与Matlab实现 (6) 5 数值解与解析解的比较 (9) 6 随时间变化的细杆上的温度分布情况 (11) 7稳定后细杆上的温度分布情况 (12) 参考文献 (13) 附录 (14)

有限长杆的一维热传导问题 问题描述 一根单位长度的细杆放入100℃的沸水中，当细杆的温度达到100℃时取出。假设细杆四周绝热；在时间t=0时，细杆两端浸入0℃的冰水中。一维热传导方 程：20t xx u a u -=，现在令21a =，从而可知本题：0t xx u u -=。现在要求细杆温度分布：(,)u x t 。 1 解析解——分离变量法 热传导偏微分方程： 0t x x u u -= (1) (0,t )(1,t )u u == (0)()u x x ?=， 其中， 001 x x ==，或 ()x ?= 100(0,1) x ∈ ， 首先令： (,)()()u x t X x T t = (2) 将(2)式带入(1)式得： ()T()()()0X x t T t X x -= 于是可得： T()()()() t X x T t X x λ==- 可以得到两个微分方程： T()()0t T t λ+= ()()0X x X x λ+= 先求解空间项： 当0λ <时， ()x x X x Ae Be λλ- --=+ 由于(0,t)(1,t)0,.u u t ==?

热传导方程及其定解问题的导出

第一章 热传导方程 本章介绍最典型的抛物型方程—热传导方程,在研究热传导,扩散等物理现象时都会遇 到这类方程. §1 热传导方程及其定解问题的导出 1.1热传导方程的导出 物理模型 在三维空间中,考虑一均匀,各向同性的物体Ω,假定它内部有热源,并且与周围介质有热交换,需要来研究物体内部温度的分布和变化. 以函数),,,(t z y x u 表示物体Ω在位置),,(z y x 及时刻t 的温度.物体内部由于各部分温度不同,产生热量的传递,它们遵循能量守恒定律. 能量守恒定律 物体内部的热量的增加等于通过物体的边界流入的热量与由物体内部的热源所生成的热量的总和 . 在物体Ω内任意截取一块D .现在时段],[21t t 上对D 使用能量守恒定律. 设),,,(t z y x u u =是温度(度),c 是比热(焦耳∕度·千克),ρ是密度(千克/米3), q 是热流密度(焦耳/秒·米2),0f 是热源强度(焦耳/千克·秒). 注意到在dt 时段内通过D 的边界D ?上小块dS 进入区域D 的热量为dSdt n q ?-(n 是 D ?的外法向),从而由能量守恒律,我们有 ,)||(21 21 120??????????+?-=-?==t t D t t D D t t t t dxdydz f dt ds n q dt dxdydz u u c ρρ （1.1） 大家知道,热量流动的原因是因为在物体内部存在温差.依据传热学中的傅立叶实验定律,在一定条件下,热流向量与温度梯度成正比 ,u k q ?-= （梯度? ?? ? ????????==?z u y u x u gradu u ,,） （1.2） 这里负号表明热量是由高温向低温流动,k 是物体的导热系数.

一维热传导方程

一维热传导方程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层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。

热传导方程

前言 本文只是针对小白而写，可以使新手对热传导理论由很浅到不浅的认识，如想更深学习热传导知识，请转其它文档。 一、概念与常量 1、温度场： 指某一时刻下，物体内各点的温度分布状态。 在直角坐标系中：； 在柱坐标系中：； 在球坐标系中：。 补充：根据温度场表达式，可分析出导热过程是几维、稳态或非稳态的现象，温度场是几维的、稳态的或非稳态的。 2、等温面与等温线： 三维物体内同一时刻所有温度相同的点的集合称为等温面； 一个平面与三维物体等温面相交所得的的曲线线条即为平面温度场中的等温线。 3、温度梯度： 在具有连续温度场的物体内，过任意一点P温度变化率最大的方向位于等温线的法线方向上。称过点P的最大温度变化率为温度梯度(temperature gradient)。用grad t表示。 定义为： 补充：温度梯度表明了温度在空间上的最大变化率及其方向，是向量，其正向与热流方向恰好相反。对于连续可导的温度场同样存在连续的温度梯度场。 在直角坐标系中： 3、导热系数 定义式：单位 导热系数在数值上等于单位温度降度(即1)下，在垂直于热流密度的单 位面积上所传导的热流量。导热系数是表征物质导热能力强弱的一个物性参数。 补充：由物质的种类、性质、温度、压力、密度以及湿度影响。 二、热量传递的三种基本方式 热量传递共有三种基本方式：热传导；热对流；热辐射 三、导热微分方程式（统一形式：） 直角坐标系：

圆柱坐标系：球坐标系：

其中，称为热扩散系数，单位，为物质密度，为物体比热容，为物体导热系数，为热源的发热率密度，为物体与外界的对流交换系数。 补充： 1处研究的对象为各向同性的、连续的、有内热源、物性参数已知的导热物体。 2稳态温度场，即则有：，此式称为泊松方程。 3无内热源的稳态温度场，则有：，此式称为拉普拉斯方程。 四、单值条件 导热问题的单值条件通常包括以下四项： 1几何条件：表示导热物体的几何形状与大小（一维、二维或三维） 2物理条件：说明导热系统的物理特性（有无内热源） 3初始条件：给出时温度场内各点温度。 数学表达式为 4边界条件：表示导热系统在边界的特征 第一类边界条件（狄利克雷边界条件）：说明物体边界的温度分布： 第二类边界条件（纽曼边界条件）：说明物体边界的热流量： 绝热边界条件 第三类边界条件（纽曼边界条件）：说明物体边界到热量与对流换热的能量平衡关系： 其中为边界处的温度，为边界的热流量，为环境温度。 五、解题步骤： 1根据具体实际问题列出导热微分方程式 2确定初始条件以及边界条件。每一维导热至少有两个边界条件，从而得到导热现象的完整数学描述。 3分析求解，得出导热物体的温度场

热传导方程的求解

应用物理软件训练 前言 MATLAB 是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。 MATLAB是矩阵实验室（Matrix Laboratory）的简称，和Mathematica、Maple 并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其

他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。 本部分主要介绍如何根据所学热传导方程的理论知识进行MATLAB数值实现可视化。

题目：热传导方程的求解 目录 一、参数说明 (1) 二、基本原理 (1) 三、MATLAB程序流程图 (3) 四、源程序 (3) 五、程序调试情况 (6) 六、仿真中遇到的问题 (9) 七、结束语 (9) 八、参考文献 (10)

一、参数说明 U=zeros(21,101) 返回一个21*101的零矩阵 x=linspace(0,1,100);将变量设成列向量 meshz(u)绘制矩阵打的三维图 axis([0 21 0 1]);横坐标从0到21，纵坐标从0到1 eps是MATLAB默认的最小浮点数精度 [X,Y]=pol2cart(R,TH);效果和上一句相同 waterfall(RR,TT,wn)瀑布图 二、基本原理 1、一维热传导问题 （1）无限长细杆的热传导定解问题 利用傅里叶变换求得问题的解是： 取得初始温度分布如下 这是在区间0到1之间的高度为1的一个矩形脉冲，于是得 （2）有限长细杆的热传导定解问题

一维非稳态导热问题的数值解

计算传热学程序报告 题目：一维非稳态导热问题的数值解 ： 学号： 学院：能源与动力工程学院 专业：工程热物理 日期：2014年5月25日

一维非稳态导热问题数值解 求解下列热传导问题： ? ?? ????=====≤≤=??- ??1,10),(,1),0(0)0,()0(01T 22ααL t L T t T x T L x t T x 1.方程离散化 对方程进行控制体积分得到： dxdt t T dxdt x T t t t e w t t t e w ? ?? ??+?+??=??α 1 2 2 ? ? -=??-???+?+e w t t t w e t t t dx T T dt x T x T )(1])()( [α 非稳态项：选取T 随x 阶梯式变化，有 x T T dx T T t p t t p e w t t t ?-=-?+?+? )()( 扩散项：选取一阶导数随时间做显示变化，有 t x T x T dt x T x T t w t e w e t t t ???-??=??-??? ?+])()[(])()[( 进一步取T 随x 呈分段线性变化，有 e P E e x T T x T )()( δ-=?? ， w W P w x T T x T )()(δ-=?? 整理可以得到总的离散方程为： 2 21x T T T t T T t W t P t E t P t t E ?+-=?-?+α 2.计算空间和时间步长 取空间步长为： h=L/N 网格Fourier 数为： 2 2 0x t x t F ??= ??= α（小于0.5时稳定）

第四章导热问题的数值解法

第四章导热问题的数值解法 1 、重点内容：①掌握导热问题数值解法的基本思路； ②利用热平衡法和泰勒级数展开法建立节点的离散方程。 2 、掌握内容：数值解法的实质。 3 、了解内容：了解非稳态导热问题的两种差分格式及其稳定性。 §4—1导热问题数值求解的基本思想及内节点方程的建立由前述 3 可知，求解导热问题实际上就是对导热微分方程在定解条件下的积分求解，从而获得分析解。但是，对于工程中几何形状及定解条件比较复杂的导热问题，从数学上目前无法得出其分析解。随着计算机技术的迅速发展，对物理问题进行离散求解的数值方法发展得十分迅速，并得到广泛应用，并形成为传热学的一个分支——计算传热学（数值传热学），这些数值解法主要有以下几种： （１）有限差分法（ 2 ）有限元方法（ 3 ）边界元方法 数值解法能解决的问题原则上是一切导热问题，特别是分析解方法无法解决的问题。如：几何形状、边界条件复杂、物性不均、多维导热问题。 一．分析解法与数值解法的异同点： ?相同点：根本目的是相同的，即确定① t=f(x ， y ， z) ；② 。 ?不同点：数值解法求解的是区域或时间空间坐标系中离散点的温度分布代替连续的温度场；分析解法求解的是连续的温度场的分布特征，而不是分散点的数值。 数值求解的基本思路及稳态导热内节点离散方程的建立 二．解法的基本概念 ?实质 对物理问题进行数值解法的基本思路可以概括为：把原来在时间、空间坐标系中连续的物理量的场，如导热物体的温度场等，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。该方法称为数值解法。 这些离散点上被求物理量值的集合称为该物理量的数值解。 2 、基本思路：数值解法的求解过程可用框图 4-1 表示。 由此可见： 1 ）物理模型简化成数学模型是基础； 2 ）建立节点离散方程是关键； 3 ）一般情况微分方程中，某一变量在某一坐标方向所需边界条件的个数等于该变量在该坐标方向最高阶导数的阶数。 ?数值求解的步骤 如图 4-2 （ a ），二维矩形域内无内热源、稳态、常物性的导热问题采用数值解法的步骤如下：（1）建立控制方程及定解条件 针对图示的导热问题，它的控制方程（即导热微分方程）为：（ a ） 边界条件： x=0 时， x=H 时， 当 y=0 时，

一维热传导方程(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 是网格界点集合。

MATLAB编辑一维热传导方程的模拟程序

求解下列热传导问题： ()()()()?????????====-=≤≤=??-??1, 10,,1,010,001222ααL t L T t T z z T L z t T z T 程序： function heat_conduction() %一维齐次热传导方程 options={'空间杆长L','空间点数N' ,'时间点数M','扩散系数alfa','稳定条件的值lambda(取值必须小于0.5)',}; topic='seting'; lines=1; def={'1','100','1000','1','0.5'}; h=inputdlg(options,topic,lines,def); L=eval(h{1}); N=eval(h{2}); M=eval(h{3}); alfa=eval(h{4}); lambda=eval(h{5});%lambda 的值必须小于0.5 %*************************************************** h=L/N;%空间步长 z=0:h:L; z=z'; tao=lambda*h^2/alfa;%时间步长 tm=M*tao;%热传导的总时间tm t=0:tao:tm; t=t'; %计算初值和边值 T=zeros(N+1,M+1); Ti=init_fun(z); To=border_funo(t); Te=border_fune(t); T(:,1)=Ti; T(1,:)=To; T(N+1,:)=Te; %用差分法求出温度T 与杆长L 、时间t 的关系 for k=1:M m=2; while m<=N T(m,k+1)=lambda*(T(m+1,k)+T(m-1,k))+(-2*lambda+1)*T(m,k); m=m+1; end;

一维热传导方程

精心整理 一维热传导方程 一． 问题介绍 考虑一维热传导方程： （1）,0),(22T t x f x u a t u ≤<+??=?? 其中a 是正常数，)(x f 是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题 1）（<∞-（2） 1）（x <0（3） (4) 二． N,M 都 三． 第 第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 τ+=-++-+-+++ )(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 利用 （8或 （9用到四．通过误差估计方程 （1）可知对任意的r ，Richardson 格式都不稳定，所以Richardson 格式绝对不稳定。 （2）当210≤r 时，向前差分格式的误差无限增长。因此向前差分格式是条件稳定。 （3）向后差分格式和六点对称格式都绝对稳定，且各自的截断误差阶分别为)(2h O +τ和)(22h O +τ。 五． 数值例子 例1令f(x)=0和a=1，可求得u （x,t ）一个解析解为u(x,t)=exp(x+t)。 1． 用向前差分格式验证得数值结果如下：

热传导方程的初值问题

§2热传导方程的初值问题 一维热传导方程的初值问题(或Cauchy 问题) ?? ???+∞<<∞-=>+∞<<∞-=??-??x x x u t x t x f x u a t u ),()0,(0 ,),,(2 2 2? （） 偏导数的多种记号xx x t u x u u x u u t u =??=??=??22,,. 问题也可记为 ?? ?+∞ <<∞-=>+∞<<∞-=-x x x u t x t x f u a u xx t ),()0,(0 ,,),(2?. Fourier 变换 我们将用Fourier 变换法求解热传导方程的柯西问题.为此我们将着重介绍Fourier 变换的基本知识.Fourier 变换在许多学科中是重要使用工具. 可积函数,设)(x f f =是定义在),(+∞-∞上的函数, 且对任意A B <，()f x 在[,]A B 上 可积，若积分 ? +∞ ∞ -dx x f )(收敛,则称)(x f 在),(+∞-∞上绝对可积。 将),(+∞-∞上绝对可积函数形成的集合记为),(1 +∞-∞L 或),(+∞-∞L ， 即{ } ∞<=+∞-∞=+∞-∞? +∞ ∞ -dx x f f L L )(| ),(),(1 ，称为可积函数空间. 连续函数空间: ),(+∞-∞上全体连续函数构成的集合,记为),(+∞-∞C , {}上连续在),(|),(+∞-∞=+∞-∞f f C ， {}上连续在),(,|),(1+∞-∞'=+∞-∞f f f C 。 定义 若),(+∞-∞∈L f ,那么积分 ),(?)(21 λπ λf dx e x f x i =? +∞ ∞ -- 有意义,称为Fourier 变换, )(? λf 称为)(x f 的Fourier 变式(或Fourier 变换的象). ? +∞ ∞ --= =dx e x f f Ff x i λπ λλ)(21)(?)( 定理 (Fourier 积分定理)若),(),(1 +∞-∞?+∞-∞∈C L f ,那么我们有

一维热传导方程的前向 、紧差分格式

中南林业科技大学 本科课程论文 学院：理学院 专业年级：09信息与计算科学一班 课程：偏微分方程数值解法 论文题目：一维热传导方程的前向Euler和紧差分格式指导教师：陈红斌 2012年7月

学生姓名：唐黎学号： 20093936分工：程序编写，数值例子 学生姓名：何雄飞学号：20093925分工：格式建立，资料收集 学生姓名：汪霄学号：20093938分工：文档编辑，资料整理 学生姓名：毛博伟学号：20093931分工：公式编辑，查找资料 学生姓名：倪新东学号：20093932分工：数据分析，查找资料 学生姓名：何凯明学号：20093924分工：数据分析，查找资料

目录 1引言 (1) 2物理背景 (1) 3网格剖分 (2) 4．1.1向前Euler格式建立 (2) 4.1.2差分格式的求解 (4) 4.1.3收敛性与稳定性 (4) 4.1.4 数值例子 (7) 4.2.1紧差分格式建立 (10) 4.2.2差分格式求解 (12) 4.2.3数值例子 (13) 总结 (17) 参考文献 (18) 附录 (19)

1 引言 本文考虑的一维非齐次热传导方程的定解问题： 22(,),0,0,u u a f x t x l t T t x ??-=<<<≤?? (,0)(),0,u x x x l φ=≤≤ (0,)(), (1,)(), 0.u t t u t t t T αβ==<≤ 其中a 为正常数，(,),(),(),()f x t x t t ?αβ为已知函数，(0)(0),(1)(0).?α?β== 目前常用的求解热传导方程的差分格式有前向Euler 差分格式、向后Euler 差分格式、Crank-Nicolson 格式、Richardson 格式[1,2,3]．本文将给出前向Euler 格式和紧差分格式，并给出其截断误差和数值例子． 2 物理背景 热传导是由于物体内部温度分布不均匀,热量要从物体内温度较高的点流向温度较低的点处.以函数(),,,u x y z t 表示物体在t 时刻,(),M M x y =处的温度,并假设 (),,u x y z 关于,,x y z 具有二阶连续偏导数,关于t 具有一阶连续偏导数.() ,,k k x y z =是物体在(),,M x y z 处的热传导系数,取正值.设物体的比热容为(),,c c x y z =,密度为 (),,x y z ρ.根据Fourier 热传导定律,热量守恒定律以及Gauss 公式得 ,u u u u c kx k k t x x y y z z ρ ????????????? =++ ? ? ???????????? ?? 如果物体是均匀的,此时,k c 以及ρ均为常数.令2 k a c ρ = ,上式方程化为 2222 2222,t u u u u a a u x y z ?? ???=++=? ?????? 若考虑物体内有热源,其热源密度函数为(),,F F x y z =,则有热源的热传导方程为 ()2,,,,t u a u f x y z t =?+ 其中F f c ρ = .