《一维导热方程-有限差分法-matlab实现23140.pdf》由会员分享，可在线阅读，更多相关《一维导热方程-有限差分法-matlab实现23140.pdf（7页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、word 专业资料-可复制编辑-欢迎下载 第五次作业（前三题写在作业纸上）一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf 文件，热扩散系数=const，22TTtx 1.用 Tylaor 展开法推导出 FTCS 格式的差分方程 2.讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。3.说明该方程的类型和定解条件，如何在程序中实现这些定解条件。4.编写 M 文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得到，添加，修改后得到。）function rechuandaopde%以下所用数据，除了 t 的范围我根据题目要求取到了 2000
2、0，其余均从 pdf 中得来 a=0.00001;%a 的取值 xspan=0 1;%x 的取值范围 tspan=0 20000;%t 的取值范围 ngrid=100 10;%分割的份数，前面的是 t 轴的，后面的是 x 轴的 f=(x)0;%初值 g1=(t)100;%边界条件一 g2=(t)100;%边界条件二 T,x,t=pdesolution(a,f,g1,g2,xspan,tspan,ngrid);%计算所调用的函数 x,t=meshgrid(x,t);mesh(x,t,T);%画图，并且把坐标轴名称改为 x，t，T xlabel(x)ylabel(t)zlabel(T)T%输出温度
3、矩阵 dt=tspan(2)/ngrid(1);%t 步长 h3000=3000/dt;h9000=9000/dt;h15000=15000/dt;%3000,9000,15000 下，温度分别在 T 矩阵的哪些行 word 专业资料-可复制编辑-欢迎下载 T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%输出三个时间下的温度分布%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面%稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x 步长 sta=4*a*dt/(dx2)*(sin(pi/2)2;if sta0
4、,sta2 fprintf(n%sn,有稳定性)else fprintf(n%sn,没有稳定性)error end%真实值计算 xe,te,Te=truesolution(a,f,g1,g2,xspan,tspan,ngrid);xe,te=meshgrid(xe,te);mesh(xe,te,Te);%画图，并且把坐标轴名称改为 xe，te，Te xlabel(xe)ylabel(te)zlabel(Te)Te%输出温度矩阵%误差计算 jmax=1/dx+1;%网格点数 rms=wuchajisuan(T,Te,jmax)rms%输出误差 function rms=wuchajisuan(T
5、,Te,jmax)for j=1:jmax word 专业资料-可复制编辑-欢迎下载 rms=(T(j)-Te(j)2/jmax)(1/2)end functionUe,xe,te=truesolution(a,f,g1,g2,xspan,tspan,ngrid)n=ngrid(1);%t 份数 m=ngrid(2);%x 份数 Ue=zeros(ngrid);xe=linspace(xspan(1),xspan(2),m);%画网格 te=linspace(tspan(1),tspan(2),n);%画网格 for j=2:n for i=2:m-1 for g=1:m-1 Ue(j,i)=
6、100-(400/(2*g-1)/pi)*sin(2*g-1)*pi*xe(j)*exp(-a*(2*g-1)2*pi2*te(i)end end end function U,x,t=pdesolution(a,f,g1,g2,xspan,tspan,ngrid)n=ngrid(1);%t 份数 m=ngrid(2);%x 份数 h=range(xspan)/(m-1);%x 网格长度 x=linspace(xspan(1),xspan(2),m);%画网格 k=range(tspan)/(n-1);%t 网格长度 t=linspace(tspan(1),tspan(2),n);%画网格 U
7、=zeros(ngrid);U(:,1)=g1(t);%边界条件 U(:,m)=g2(t);U(1,:)=f(x);%初值条件%差分计算 for j=2:n for i=2:m-1 word 专业资料-可复制编辑-欢迎下载 U(j,i)=(1-2*a*k/h2)*U(j-1,i)+a*k/h2*U(j-1,i-1)+a*k/h2*U(j-1,i+1);end end 5.将温度随时间变化情况用曲线表示 00.20.40.60.8100.511.52x 104020406080100 xtT 6.给出 3000、9000、15000 三个时刻的温度分布情况，对温度随时间变化规律做说明。T3000
8、=100.0000 63.4362 34.2299 15.8021 7.4641 7.4641 15.8021 34.2299 63.4362 100.0000 T9000=100.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.6839 65.6076 81.6930 100.0000 T15000=100.0000 89.9415 81.0962 74.5310 71.0378 71.0378 74.5310 81.0962 89.9415 100.0000 根据数据分析，在同一个 x 点上温度随时间的增加而增加，但增幅变小。x-T 图形仍为
9、抛物线，但随着时间的增加，极值变小，图像变得平缓。7.用计算数据说明，并结合差分方程余项，空间、时间间隔对求解精度影响。数据量较大，且原理相同，我取一个向量演示一下。保持空间间隔不变，修改时间间隔，时间间隔加大，得到的误差加大。保持时间间隔不变，修改空间间隔，空间间隔加大，得到的误差加大。修改空间间隔的误差在增量比修改时间间隔的大。word 专业资料-可复制编辑-欢迎下载 从方差余项上来看，（没有公式编辑器。只能从 ppt 里粘贴了）这个余项里的t，x 都在分母上，所以与误差成正比，且x 的次数应该是比t 高，故影响较大。8.用计算数据说明，稳定性要求对求解精度的影响。修改稳定性，即修改 x
10、和 t 分的份数（ngrid），之后看误差。稳定性越高，解的精度越高。即在满足稳定性要求（a*t/(x2)0.5）时，a*t/(x2)越接近 0，误差越小。从概念上理解，稳定性越好，对引入时间层误差的抑制能力越强。所以误差越小。二、调用 MATLAB 函数完成上述计算 1.编写 M 文件求解上述方程，并用适当的文字对程序做出说明。function pdepediaoyong m=0;x=linspace(0,1,11);%x 的网格 t=linspace(0,20000,101);%t 的网格 sol=pdepe(m,pdefun,icfun,bcfun,x,t);%调用函数 T=sol(:,
11、:,1);%解 figure;%画图 surf(x,t,T)xlabel(x)ylabel(t)zlabel(T)dt=20000/100;%t 步长 h3000=3000/dt;h9000=9000/dt;h15000=15000/dt;%3000,9000,15000 下，温度分别在 T 矩阵的哪些行 T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%输出三个时间下的温度分布%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面 word 专业资料-可复制编辑-欢迎下载 function c,f,s=pdefun(x,t,T,Du
12、Dx)%PDE 方程函数 c=100000;f=DuDx;s=0;function u0=icfun(x)%初始条件函数 u0=0;function pl,ql,pr,qr=bcfun(xl,Tl,xr,Tr,t)%边界条件函数 pl=Tl-100;ql=0;pr=Tr-100;qr=0;2.将温度随时间变化情况用曲线表示。00.20.40.60.8100.511.52x 104020406080100120 3.给出 3000、9000、15000 三个时刻的温度分布情况。T3000=100.0000 67.1058 39.8611 21.1973 10.9885 7.8279 10.988
13、5 21.1973 39.8611 67.1058 100.0000 word 专业资料-可复制编辑-欢迎下载 T9000=100.0000 83.4839 68.6032 56.8191 49.2705 46.6732 49.2705 56.8191 68.6032 83.4839 100.0000 T15000=100.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 函数的稳定性要用什么检验。傅里叶级数检验不适用。