一维热传导方程

1、精心整理一维热传导方程一 问题介绍考虑一维热传导方程：（1）其中a是正常数，是给定的连续函数。按照定解条件的不同给法，可将方程（1）的定解问题分为两类：第一类、初值问题（也称Cauthy问题）：求具有所需次数偏微商的函数，满足方程（1）（）和初始条件：（2）第二类、初边值问题（也称混合问题）：求具有所需次数偏微商的函数，满足方程（1）（）和初始条件：（3）及边值条件(4)假定在相应区域光滑，并且在满足相容条件，使上述问题有唯一充分光滑的解。二 区域剖分考虑边值问题（1），（4）的差分逼近。去空间步长和时间步长，其中N,M都是正整数。用两族平行直线：将矩形域分割成矩形网格，网格节点为。以表示网格

2、内点集合，即位于开矩形G的网点集合；表示所有位于闭矩形的网点集合；=-是网格界点集合。三 离散格式第k+1层值通过第k层值明显表示出来，无需求解线性代数方程组，这样的格式称为显格式。第k+1层值不能通过第k层值明显表示出来，而由线性代数方程组确定，这样的格式称为隐格式。1. 向前差分格式（5），其中j=1,2,，N-1，k=1,2,M-1。以表示网比。则方程（5）可以改写为：易知向前差分格式是显格式。2. 向后差分格式（6），其中j=1,2,，N-1，k=1,2,M-1，易知向前差分格式是显格式。3. 六点对称格式（Grank-Nicolson格式）将向前差分格式和向后差分格式作算术平均，即得

3、到六点对称格式：（7）=利用和边值便可逐层求到。六点对称格式是隐格式，由第k层计算第k+1层时需解线性代数方程组（因系数矩阵严格对角占优，方程组可唯一求解）。将其截断误差于（=）展开，则得=。4. Richardson格式（8），或（9）。这是三层显示差分格式。截断误差阶为。为了使计算能够逐层进行，除初值外，还要用到，这可以用前述二层差分格式计算（为保证精度，可将0,分成若干等份）。四 格式稳定性通过误差估计方程（1）可知对任意的r，Richardson格式都不稳定，所以Richardson格式绝对不稳定。（2）当时，向前差分格式趋于稳定；当时，向前差分格式的误差无限增长。因此向前差分格式是条

4、件稳定。（3）向后差分格式和六点对称格式都绝对稳定，且各自的截断误差阶分别为和。五 数值例子例1令f(x)=0和a=1，可求得u（x,t）一个解析解为u(x,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.11

5、11112.1766302.1737190.0029110.3333330.1666671.6487211.6456190.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.943620

6、0.0041140.6666670.3333332.7182822.7136510.0046300.3333330.3888892.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.55

7、55562.4324252.4272330.0051930.6666670.5555563.3947233.3888780.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.003483

8、0.0069080.3333330.7777783.0377323.0312430.0064880.6666670.7777784.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.94

9、44445.0083734.9997450.008628当n等于2和m等于17时最大误差为0.008628其中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.356

10、4182.3653460.0089270.8571432.7182822.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层

11、结果时间节点Tk=0.4285710.1428571.7707951.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.3740

12、620.0176440.4285712.7182822.7404570.0221750.5714293.1357153.1590580.0233430.7142863.6172513.6378270.0205760.8571434.1727344.1858510.013117第4层的最大误差是0.023343第5层结果时间节点Tk=0.7142860.1428572.3564182.3682760.0118570.2857142.7182822.7388800.0205980.4285713.1357153.1616000.0258860.5714293.6172513.6444850.027

13、2340.7142864.1727344.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层的

14、最大误差是0.031544当n等于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.7182822.

15、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.7707951.77

16、12330.0004380.2857142.0427272.0434760.0007490.4285712.3564182.3573550.0009370.5714292.7182822.7192680.0009860.7142863.1357153.1365920.0008770.8571433.6172513.6178250.000574第3层的最大误差是0.000986第4层结果时间节点Tk=0.5714290.1428572.0427272.0432220.0004950.2857142.3564182.3572870.0008690.4285712.7182822.7193770.0

17、010950.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.1727344.1739000.0011660.

18、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当n等于6时最大误差为0.0015314用Ric

19、hardson格式验证数值结果如下：请输入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.7182822.7189490.0006670.1666670.3333331.6487211

20、.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.0404110.3333330.5000002.3009762.2916210.0093550.5000000.5000002.7182822.7075140.0107680.666

21、6670.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.4816891.5458782.9358110.1666670.8333332.7182823

22、5.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.h>#defineMax_N1000doubleuMax_NMax_N,bMax_N,

23、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);un1=yn;for(i=n-1;i>0;i-)ui1=yi-betai*ui+11;intmain()/

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

25、(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+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

26、=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("%d",&n)printf("请输入m的值(输入0结束程序):");system("PASUE");return0;附录2

27、#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMax_N1000doubleuMax_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<

28、;=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");if(scanf("%d",&n)printf("请输入m的值(输入0结束程序):n");while(scanf("%d",&m)&&a

29、mp;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;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+

30、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);elsefi=ui;catchup();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=te

31、m+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);/printf("当n等于%d时其二范数为%lfn",n,Errj);/if(j!=0)printf("log(Errj-1/Errj)=%lfnn",log(Errj-1/Errj);/j+;printf("请输入n的值(输入0结束程序

32、):");if(scanf("%d",&n)printf("请输入m的值(输入0结束程序):");system("PASUE");return0;附录3#include<stdio.h>#include<stdlib.h>#include<math.h>#defineMax_N1000doubleuMax_N,bMax_N,aMax_N,cMax_N,fMax_N,errMax_N,xMax_N,yMax_N,betaMax_N,ErrMax_N;intn,m;/将空间区间【0，1

33、】分为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);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&qu

34、ot;);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;for(k=1;k<=m

35、;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)*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

36、/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);elsefi=r/2*ui+1+(1-r)*ui+r/2*ui-1;catchup();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

37、)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);/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_N1000doubleuMax_NMax_N,bMax_N,aMax_N,cMax_N,fMax_N,errMax_NMax_N,xMax_NMax_N,yM