ftcs有限差分法c语言程序,一维导热方程-有限差分法-matlab实现.docx

一维导热方程-有限差分法-matlab实现.docx

第五次作业(前三题写在作业纸上)一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf文件，热扩散系数const，1. 用Tylaor展开法推导出FTCS格式的差分方程2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。4. 编写M文件求解上述方程，并用适当的文字对程序做出说明。(部分由网络搜索得到，添加，修改后得到。)function rechuandaopde以下所用数据，除了t的范围我根据题目要求取到了20000，其余均从pdf中得来a0.00001;a的取值xspan0 1;x的取值范围tspan0 20000;t的取值范围ngrid100 10;分割的份数，前面的是t轴的，后面的是x轴的fx0;初值g1t100;边界条件一g2t100;边界条件二T,x,tpdesolutiona,f,g1,g2,xspan,tspan,ngrid;计算所调用的函数x,tmeshgridx,t;meshx,t,T;画图，并且把坐标轴名称改为x，t，TxlabelxylabeltzlabelTT输出温度矩阵dttspan2/ngrid1;t步长h30003000/dt;h90009000/dt;h1500015000/dt;3000,9000,15000下，温度分别在T矩阵的哪些行T3000Th3000,T9000Th9000,T15000Th15000,输出三个时间下的温度分布不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面稳定性讨论,傅里叶级数法dxxspan2/ngrid2;x步长sta4*a*dt/dx2*sinpi/22;if sta0,sta2 fprintfnsn,有稳定性else fprintfnsn,没有稳定性 errorend真实值计算xe,te,Tetruesolutiona,f,g1,g2,xspan,tspan,ngrid;xe,temeshgridxe,te;meshxe,te,Te;画图，并且把坐标轴名称改为xe，te，TexlabelxeylabeltezlabelTeTe输出温度矩阵误差计算jmax1/dx1;网格点数rmswuchajisuanT,Te,jmaxrms输出误差function rmswuchajisuanT,Te,jmaxfor j1jmax rmsTj-Tej2/jmax1/2endfunctionUe,xe,tetruesolutiona,f,g1,g2,xspan,tspan,ngridnngrid1;t份数mngrid2;x份数Uezerosngrid;xelinspacexspan1,xspan2,m;画网格 telinspacetspan1,tspan2,n;画网格for j2n for i2m-1 for g1m-1Uej,i100-400/2*g-1/pi*sin2*g-1*pi*xej*exp-a*2*g-12*pi2*tei end endendfunction U,x,tpdesolutiona,f,g1,g2,xspan,tspan,ngridnngrid1;t份数mngrid2;x份数hrangexspan/m-1;x网格长度 xlinspacexspan1,xspan2,m;画网格 krangetspan/n-1; t网格长度tlinspacetspan1,tspan2,n;画网格Uzerosngrid;U,1g1t;边界条件U,mg2t;U1,fx;初值条件 差分计算 for j2n for i2m-1 Uj,i1-2*a*k/h2*Uj-1,ia*k/h2*Uj-1,i-1a*k/h2*Uj-1,i1; endend5. 将温度随时间变化情况用曲线表示6. 给出3000、9000、15000三个时刻的温度分布情况，对温度随时间变化规律做说明。T3000100.0000 63.4362 34.2299 15.8021 7.4641 7.4641 15.8021 34.2299 63.4362 100.0000T9000100.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.6839 65.6076 81.6930 100.0000T15000100.0000 89.9415 81.0962 74.5310 71.0378 71.0378 74.5310 81.0962 89.9415 100.0000根据数据分析，在同一个x点上温度随时间的增加而增加，但增幅变小。x-T图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。7. 用计算数据说明，并结合差分方程余项，空间、时间间隔对求解精度影响。数据量较大，且原理相同，我取一个向量演示一下。保持空间间隔不变，修改时间间隔，时间间隔加大，得到的误差加大。保持时间间隔不变，修改空间间隔，空间间隔加大，得到的误差加大。修改空间间隔的误差在增量比修改时间间隔的大。从方差余项上来看，(没有公式编辑器只能从ppt里粘贴了)这个余项里的t，x都在分母上，所以与误差成正比，且x的次数应该是比t高，故影响较大。8. 用计算数据说明，稳定性要求对求解精度的影响。修改稳定性，即修改x和t分的份数(ngrid)，之后看误差。稳定性越高，解的精度越高。即在满足稳定性要求(a*t/x20.5)时，a*t/x2越接近0，误差越小。从概念上理解，稳定性越好，对引入时间层误差的抑制能力越强。所以误差越小。二、调用MATLAB函数完成上述计算1. 编写M文件求解上述方程，并用适当的文字对程序做出说明。function pdepediaoyongm0;xlinspace0,1,11;x的网格tlinspace0,20000,101;t的网格sol pdepem,pdefun,icfun,bcfun,x,t;调用函数Tsol,,1;解figure;画图surfx,t,TxlabelxylabeltzlabelTdt20000/100;t步长h30003000/dt;h90009000/dt;h1500015000/dt;3000,9000,15000下，温度分别在T矩阵的哪些行T3000Th3000,T9000Th9000,T15000Th15000,输出三个时间下的温度分布不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面function c,f,spdefunx,t,T,DuDxPDE方程函数c100000;fDuDx;s0;function u0icfunx初始条件函数u00;function pl,ql,pr,qrbcfunxl,Tl,xr,Tr,t边界条件函数pl Tl-100;ql 0;pr Tr-100;qr 0;2. 将温度随时间变化情况用曲线表示。3. 给出3000、9000、15000三个时刻的温度分布情况。T3000 100.0000 67.1058 39.8611 21.1973 10.9885 7.8279 10.9885 21.1973 39.8611 67.1058 100.0000T9000100.0000 83.4839 68.6032 56.8191 49.2705 46.6732 49.2705 56.8191 68.6032 83.4839 100.0000T15000100.0000 90.8310 82.5601 75.9972 71.7845 70.3330 71.7845 75.9972 82.5601 90.8310 100.0000根据数据分析，在同一个x点上温度随时间的增加而增加，但增幅变小。x-T图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。4. 用计算数据说明，空间、时间间隔对求解精度影响，并与有限差分法的计算结果做比较。调用前面做出来的真实值,跟pdepe做出来的值计算误差，再与有限差分法的误差比较。用pdepe函数求的误差远小于有限差分法，所以pdepe函数法更精确。5. 用计算数据说明，有无稳定性要求，为什么若有，如何对求解精度的影响。不知道这个pdepe函数的稳定性要用什么检验。傅里叶级数检验不适用。