一维热传导方程

1、一 维 热 传 导 方 程1. 问题介绍考虑一维热传导方程：2a2 f(x),0 t T, t x其中a是正常数，f (x)是给定的连续函数。按照定解条件的不同给法，可将方程(1)的定解问题分为两类：第一类、初值问题 (也称Cauthy问题)：求具有所需次数偏微商的函数u(x,t),满足方程(1)( x )和初始条件：(2) U(x,0)(x), x第二类、初边值问题(也称混合问题)：求具有所需次数偏微商的函数u(x,t),满足方程(1) (0 x l )和初始条件：(3) U(x,0)(x),0 x l及边值条件u(0,t) u(l,t) 0. 0 t T假定(x)在相应区域光滑，并且在 x

2、 0,l满足相容条件，使上述问题有唯一充分光滑的解。2. 区域剖分考虑边值问题(1) , (4)的差分逼近。去空间步长 h l/N和时间步长T/M ,其中N,M都是正整数。用两族平行直线：将矩形域 G 0 x l;0 t T分割成矩形网格，网格节点为(xj,tk)。以Gh表示网格内点集合，即位于开矩形 G的网点集合；Gh表示所有位于闭矩形 G的网点集合；h = Gh-Gh是网格界点集合。3. 离散格式第k+1层值通过第k层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格 式。第k+1层值不能通过第 k层值明显表示出来，而由线性代数方程组确定，这样的格式称为 隐格式。1 .向前差分格式

3、k 1 kkk ku j u j u j 1 2u j u i 1-L a-j fj,h2jfjf(Xj),u0 j(Xj), ukuN0,其中 j=1,2,，N-1 , k=1,2，,M-1。以 r2 .a /h表示网比。则方程(5)可以改写为:易知向前差分格式是显格式2 .向后差分格式(6)k 1k 1 k 1 krUj i (1 2u)Uj ru j i Uj fj,.0k . ku j j(Xj )，u0uN0 ,其中j=1,2，N-1, k=1,2，,M-1 ,易知向前差分格式是显格式。3 . 六点对称格式(Grank-Nicolson格式)将向前差分格式和向后差分格式作算术平均，即

4、得到六点对称格式:juki (1 r)uk1r kuj2 jr k= -uj2 jk r k(1r)uj2uj 1fj利用u；和边值便可逐层求到kuj。六点对称格式是隐格式，由第 k层计算第k+1层时需解线性代数方程组(因系数矩阵严格对角占优，方程组可唯一求解)将其截断误差R：(x)于 区立1) jk .， 八 1、(t 1 = (k -)展开，则得k22lRk(x) = O( 2 h2)4 . Richardson 格式(8)k 1 k 1uj ujkuja -2u：h2kuj 1fj，k 1 k k k k 1uj2r(uj 1 2uj uj 1) uj2 fj，22、 0这是三层显示差分

5、格式。截断误差阶为0( 2 h2) 0为了使计算能够逐层进行，除初值u；1外，还要用到uj ,这可以用前述二层差分格式计算(为保证精度，可将0,分成若干等份)4. 格式稳定性通过误差估计方程11)可知对任意的 r, Richardson格式都不稳定，所以 Richardson格式绝对不稳定。1 1(2)当0 r 一时，向前差分格式趋于稳定；当 r 一时，向前差分格式的误差无限增长。22因此向前差分格式是条件稳定。3 ）向后差分格式和六点对称格式都绝对稳定，且各自的截断误差阶分别为O( h2)和O( 2 h2)。五 数值例子例1令f(x)=0和a=1 ,可求得u (x,t) 一个解析解为 u(x

6、,t)=exp(x+t)1 用向前差分格式验证得数值结果如下：请输入 n 的值 (输入 0 结束程序):2请输入 m 的值 ( 输入 0 结束程序):17xjtk 真实值 xik 近似值 uik 误差 errik0.3333330.0555561.4753411.4738670.0014740.6666670.0555562.0590042.0569470.0020570.3333330.1111111.5596231.5570370.0025860.6666670.1111112.1766302.1737190.0029110.3333330.1666671.6487211.6456190.

7、0031020.6666670.1666672.3009762.2973850.0035910.3333330.2222221.7429091.7393730.0035360.6666670.2222222.4324252.4284450.0039810.3333330.2777781.8424771.8386470.0038310.6666670.2777782.5713842.5670480.0043370.3333330.3333331.9477341.9436200.0041140.6666670.3333332.7182822.7136510.0046300.3333330.3888

8、892.0590042.0546320.0043720.6666670.3888892.8735712.8686440.0049270.3333330.4444442.1766302.1719920.0046380.6666670.4444443.0377323.0325120.0052200.3333330.5000002.3009762.2960680.0049080.6666670.5000003.2112713.2057440.0055260.3333330.5555562.4324252.4272330.0051930.6666670.5555563.3947233.3888780.

9、0058450.3333330.6111112.5713842.5658940.0054910.6666670.6111113.5886563.5824750.0061810.3333330.6666672.7182822.7124760.0058050.6666670.6666673.7936683.7871330.0065350.3333330.7222222.8735712.8674340.0061370.6666670.7222224.0103924.0034830.0069080.3333330.7777783.0377323.0312430.0064880.6666670.7777

10、784.2394964.2321930.0073030.3333330.8333333.2112713.2044110.0068590.6666670.8333334.4816894.4739690.0077210.3333330.8888893.3947233.3874720.0072510.6666670.8888894.7377184.7295560.0081620.3333330.9444443.5886563.5809910.0076650.6666670.9444445.0083734.9997450.008628当 n 等于 2 和 m 等于 17 时最大误差为0.008628其

11、中 r=1/2 ，格式是稳定的。2 .用向后差分格式验证得数值结果如下：请输入 n 的值（输入0 结束程序）:6请输入 m 的值 （ 输入 0 结束程序）:6xj 真实值 xi 近似值 ui 误差 erri第 1 层结果时间节点Tk=0.1428570.1428571.3307121.3350020.0042890.2857141.5350631.5423580.0072950.4285711.7707951.7799490.0091540.5714292.0427272.0525240.0097970.7142862.3564182.3653460.0089270.8571432.71828

12、22.7242560.005974第1 层的最大误差是0.009797第 2 层结果时间节点Tk=0.2857140.1428571.5350631.5417970.0067340.2857141.7707951.7824240.0116290.4285712.0427272.0573450.0146180.5714292.3564182.3718950.0154770.7142862.7182822.7320700.0137880.8571433.1357153.1446330.008919第 2 层的最大误差是0.015477第3 层结果时间节点Tk=0.4285710.1428571.7

13、707951.7793240.0085290.2857142.0427272.0575180.0147910.4285712.3564182.3750100.0185920.5714292.7182822.7378840.0196020.7142863.1357153.1530410.0173260.8571433.6172513.6283370.011086第3 层的最大误差是0.019602第 4 层结果时间节点Tk=0.5714290.1428572.0427272.0528880.0101610.2857142.3564182.3740620.0176440.4285712.71828

14、22.7404570.0221750.5714293.1357153.1590580.0233430.7142863.6172513.6378270.0205760.8571434.1727344.1858510.013117第 4 层的最大误差是0.023343第5 层结果时间节点Tk=0.714286 0.1428572.3564182.3682760.0118570.2857142.7182822.7388800.0205980.4285713.1357153.1616000.0258860.5714293.6172513.6444850.0272340.7142864.1727344.

15、1967160.0239820.8571434.8135204.8287880.015268第5 层的最大误差是0.027234第6 层结果时间节点Tk=0.8571430.1428572.7182822.7320170.0137350.2857143.1357153.1595780.0238640.4285713.6172513.6472400.0299890.5714294.1727344.2042780.0315440.7142864.8135204.8412870.0277670.8571435.5527085.5703780.017670第 6 层的最大误差是0.031544当 n

16、等于 6 时最大误差为0.0315443 .用六点对称格式验证数值结果如下：请输入 n 的值（输入0 结束程序）:6请输入 m 的值 （ 输入 0 结束程序）:6xj 真实值 xi 近似值 ui 误差 erri第 1 层结果时间节点Tk=0.1428570.1428571.3307121.3309880.0002760.2857141.5350631.5355220.0004590.4285711.7707951.7713680.0005730.5714292.0427272.0433500.0006230.7142862.3564182.3570040.0005850.8571432.718

17、2822.7186920.000410第1 层的最大误差是0.000623第 2 层结果时间节点Tk=0.2857140.1428571.5350631.5354240.0003610.2857141.7707951.7714350.0006400.4285712.0427272.0435380.0008100.5714292.3564182.3572680.0008490.7142862.7182822.7190160.0007340.8571433.1357153.1361620.000447第 2 层的最大误差是0.000849第3 层结果时间节点Tk=0.4285710.1428571

18、.7707951.7712330.0004380.2857142.0427272.0434760.0007490.4285712.3564182.3573550.0009370.5714292.7182822.7192680.0009860.7142863.1357153.1365920.000877 0.8571433.6172513.6178250.000574第3 层的最大误差是0.000986第 4 层结果时间节点Tk=0.5714290.1428572.0427272.0432220.0004950.2857142.3564182.3572870.0008690.4285712.71

19、82822.7193770.0010950.5714293.1357153.1368670.0011530.7142863.6172513.6182610.0010100.8571434.1727344.1733650.000631第 4 层的最大误差是0.001153第5 层结果时间节点Tk=0.7142860.1428572.3564182.3569990.0005810.2857142.7182822.7192840.0010030.4285713.1357153.1369730.0012580.5714293.6172513.6185730.0013220.7142864.172734

20、4.1739000.0011660.8571434.8135204.8142710.000751第5 层的最大误差是0.001322第6 层结果时间节点Tk=0.8571430.1428572.7182822.7189450.0006630.2857143.1357153.1368710.0011570.4285713.6172513.6187050.0014550.5714294.1727344.1742650.0015310.7142864.8135204.8148670.0013470.8571435.5527085.5535590.000851第 6 层的最大误差是0.001531当

21、n 等于 6 时最大误差为0.0015314用Richardson 格式验证数值结果如下：请输入 n 的值 （输入 0 结束程序）:5请输入 m 的值 （ 输入 0 结束程序）:5xjtk 真实值 xik 近似值 uik 误差 errik0.1666670.1666671.3956121.3960800.0004680.3333330.1666671.6487211.6494810.0007600.5000000.1666671.9477341.9486490.0009150.6666670.1666672.3009762.3018880.0009120.8333330.1666672.718

22、2822.7189490.0006670.1666670.3333331.6487211.6455400.0031820.3333330.3333331.9477341.9448060.0029280.5000000.3333332.3009762.2975830.0033930.6666670.3333332.7182822.7136000.0046820.8333330.3333333.2112713.2040990.0071710.1666670.5000001.9477341.9881450.040411 0.3333330.5000002.3009762.2916210.009355

23、0.5000000.5000002.7182822.7075140.0107680.6666670.5000003.2112713.1956850.0155860.8333330.5000003.7936683.9077790.1141110.1666670.6666672.3009761.2141561.0868200.3333330.6666672.7182823.2938220.5755410.5000000.6666673.2112713.1649070.0463640.6666670.6666673.7936685.4006921.6070240.8333330.6666674.48

24、16891.5458782.9358110.1666670.8333332.71828235.74707433.0287920.3333330.8333333.211271-24.21136127.4226310.5000000.8333333.79366831.08392727.2902590.6666670.8333334.481689-69.89150974.3731980.8333330.8333335.29449095.14889189.854401附录 1#include<stdio.h>#include<stdlib.h>#include<math.

25、h>#defineMax_N1000double uMax_NMax_N,bMax_N,aMax_N,cMax_N,fMax_N, errMax_NMax_N,xMax_NMax_N,yMax_N,betaMax_N,ErrMax_N;intn,m;/ 将空间区间【0， 1 】分为 n 等份；时间区间【0， 1 】分为 m 等份voidcatchup() inti;beta1=c1/b1;for(i=2;i<n;i+)betai=ci/(bi-ai*betai-1);y1=f1/b1;for(i=2;i<=n;i+)yi=(fi-ai*yi-1)/(bi-ai*betai-1

26、);un1=yn;for(i=n-1;i>0;i-)ui1=yi-betai*ui+11; intmain()/ 一维热传导方程的向前差分格式intk,i;doubleh,t,r;doublepi=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)

27、,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+)/ 初值条件ui0=exp(i*h);for(k=0;k<=m+1;k+)/ 边值条件u0k=exp(k*t);un+1k=exp(n+1)*h+k*t);printf("xjtk 真实值 xik 近似值 uik 误差 errikn");for(k=1;k<=m;k+)for(intj=1;j<=n;j+)ujk=r*uj+1k-1+(1-2*r)*ujk-1

28、+r*uj-1k-1;doubleT=0;for(k=1;k<=m;k+)for(i=1;i<=n;i+)xik=exp(i*h+k*t);errik=xik>uik?xik-uik:uik-xik;printf("%lf%lf%lf%lf%lfn",i*h,k*t,xik,uik,errik); if(errik>T)T=errik;printf(" 当 n 等于 %d 和 m 等于 %d 时最大误差为%lfn",n,m,T);printf(" 请输入 n 的值(输入0 结束程序):");if(scanf(&

29、quot;%d",&n)printf(" 请输入 m 的值(输入0 结束程序):");system("PASUE");return0;附录 2#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMax_N1000double uMax_N,bMax_N,aMax_N,cMax_N,fMax_N,errMax_N,xMax_N,yMax_N,betaMax_N,ErrMax_N;intn,m;/ 将空间区间【0，1 】分为 n 等份；时间区

30、间【0， 1 】分为 m 等份voidcatchup()inti;beta1=c1/b1;for(i=2;i<n;i+)betai=ci/(bi-ai*betai-1);y1=f1/b1;for(i=2;i<=n;i+)yi=(fi-ai*yi-1)/(bi-ai*betai-1);un=yn;for(i=n-1;i>0;i-)ui=yi-betai*ui+1;intmain()/ 一维热传导方程的六点对称格式intk,i;doubleh,t,r;/intj=0;doublepi=3.1415627;printf(" 请输入 n 的值 (输入 0 结束程序):n&q

31、uot;);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 真实值 xi 近似值 ui 误差 errin");doubleM=0;/doubletemp=0;

32、for(k=1;k<=m;k+) b1=1+2*r;c1=-r;an=-r;bn=1+2*r;if(k=1)f1=exp(h)+r*exp(t);fn=exp(n*h)+r*exp(n+1)*h+t);elsef1=u1+r*exp(k*t);fn=un+r*exp(n+1)*h+k*t);for(i=2;i<n;i+)bi=1+2*r;ai=-r;ci=-r;if(k=1) fi=exp(i*h);else fi=ui;catchup();printf(" 第 %d 层结果时间节点Tk=%lfn",k,k*t);doubleT=0;/doubletem=0;f

33、or(i=1;i<=n;i+)xi=exp(i*h+k*t);erri=xi>ui?xi-ui:ui-xi;printf("%lf%lf%lf%lfn",i*h,xi,ui,erri);if(erri>T)T=erri;/tem=tem+erri*erri*h;printf(" 第 %d 层的最大误差是%lfn",k,T);if(T>M)M=T;/temp=temp+tem;printf("n");/Errj=temp;printf(" 当 n 等于 %d 时最大误差为%lfn",n,M);

34、/printf(" 当 n 等于 %d 时其二范数为%lfn",n,Errj);/if(j!=0)printf("log(Errj-1/Errj)=%lfnn",log(Errj-1/Errj);/j+;printf(" 请输入 n 的值(输入0 结束程序):");if(scanf("%d",&n)printf(" 请输入 m 的值(输入0 结束程序):");system("PASUE");return0;附录 3#include<stdio.h> #in

35、clude<stdlib.h>#include<math.h>#defineMax_N1000double uMax_N,bMax_N,aMax_N,cMax_N,fMax_N,errMax_N, xMax_N,yMax_N,betaMax_N,ErrMax_N;intn,m;/ 将空间区间【0，1 】分为 n 等份；时间区间【0， 1 】分为 m 等份voidcatchup()inti;beta1=c1/b1;for(i=2;i<n;i+)betai=ci/(bi-ai*betai-1);y1=f1/b1;for(i=2;i<=n;i+)yi=(fi-ai

36、*yi-1)/(bi-ai*betai-1);un=yn;for(i=n-1;i>0;i-)ui=yi-betai*ui+1;intmain()/ 一维热传导方程的六点对称格式intk,i;doubleh,t,r;/intj=0;doublepi=3.1415627;printf(" 请输入 n 的值 (输入 0 结束程序):n");if(scanf("%d",&n)printf(" 请输入 m 的值(输入0 结束程序):n");while(scanf("%d",&m)&&m&

37、amp;&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 真实值 xi 近似值 ui 误差 errin");doubleM=0;/doubletemp=0;for(k=1;k<=m;k+) b1=1+r;c1=-r/2;an=-r/2;bn=1+r;if(k=1) f1=r/2*exp(2*h)+(1-r)*exp(h)+r/2+r/2*exp(t);fn=r/2*exp(n+1)*h)+(1-r

38、)*exp(n*h)+r/2*exp(n-1)*h)+r/2*exp(n+1)*h+t);elsef1=r/2*u2+(1-r)*u1+r/2*exp(k-1)*t)+r/2*exp(k*t);fn=r/2*un-1+(1-r)*un+r/2*exp(n+1)*h+(k-1)*t)+r/2*exp(n+1)*h+k*t); for(i=2;i<n;i+) bi=1+r;ai=-r/2;ci=-r/2;if(k=1) fi=r/2*exp(i+1)*h)+(1-r)*exp(i*h)+r/2*exp(i-1)*h);else fi=r/2*ui+1+(1-r)*ui+r/2*ui-1; c

39、atchup();printf(" 第 %d 层结果时间节点Tk=%lfn",k,k*t);doubleT=0;doubletem=0;for(i=1;i<=n;i+) xi=exp(i*h+k*t);erri=xi>ui?xi-ui:ui-xi;printf("%lf%lf%lf%lfn",i*h,xi,ui,erri);if(erri>T)T=erri;/tem=tem+erri*erri*h; printf(" 第 %d 层的最大误差是%lfn",k,T);if(T>M)M=T;/temp=temp+te

40、m;printf("n");/Errj=temp;printf(" 当 n 等于 %d 时最大误差为%lfn",n,M);/printf(" 当 n 等于 %d 时其二范数为%lfn",n,Errj);/if(j!=0)printf("log(Errj-1/Errj)=%lfnn",log(Errj-1/Errj);/j+;printf(" 请输入 n 的值(输入0 结束程序):");if(scanf("%d",&n)printf(" 请输入 m 的值 (输入 0 结束程序):");system("PASUE");return0;附录 4#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMax_N1000double uMax_NMax_N,bMax_N,aMax_N,cMax_