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

1、第五次作业(前三题写在作业纸上)一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf文件,热扩散系数a =const,:T；：2T=2:t: x1. 用Tylaor展开法推导出FTCS式的差分方程2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。4. 编写M文件求解上述方程，并用适当的文字对程序做出说明。(部分由网络搜索得至L添加，修改后得到。)function rechuandaopde%以下所用数据，除了t的范围我根据题目要求取到了20000，其余均从pdf中得来a=0.00001

2、;%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, Txlabel('x')ylabel('t')zlabel('T')T%输出温度矩阵d

3、t=tspan(2)/ngrid(1);%t 步长h3000=3000/dt;h9000=9000/dt;h15000=15000/dt;%3000,9000,15000下，温度分别在 T矩阵的哪些行T3000=T(h3000,:)T9000=T(h9000,:)T15000=T(h15000,:)%输出三个时间下的温度分布%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面%稳定性讨论，傅里叶级数法dx=xspan(2)/ngrid(2);%x 步长sta=4*a*dt/(dxA2)*(sin(pi/ 2)A2;if sta>0,sta<2fprintf('n%

4、sn','有稳定性')elsefprintf('n%sn','没有稳定性')errorend%真实值计算xe,te,Te=truesolution(a,f,g1,g2,xspan,tspan,ngrid);xe,te=meshgrid(xe,te);mesh(xe,te,Te);%画图，并且把坐标轴名称改为xe, te, Texlabel('xe')ylabel('te')zlabel('Te')Te%谕出温度矩阵%误差计算jmax=1/dx+1;% 网格点数rms=wuchajisuan(

5、T,Te,jmax)rms%输出误差function rms=wuchajisuan(T,Te,jmax)for j=1:jmaxrms=(T(j)-Te(j)A2/jmax)A(1 /2)endfunctionUe,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:nfor i=2:m-1

6、for g=1:m-1Ue(j,i)=100-(400/(2*g-1)/pi)*sin(2*g-1)*pi*xe(j)*exp(-a*(2*g-1)A2*piA2*te(i)endendendfunction 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),tsp

7、an,n);% 画网格U=zeros(ngrid);U(:,1)=g1(t);% 边界条件U(:,m)=g2(t);U(1,:)=f(x);% 初值条件%差分计算for j=2:nfor i=2:m-1U(j,i)=(1-2*a*k/hA2)*U(j-1,i)+a*k/hA2*U(j-1,i-1)+a*k/hA2*U(j-1,i+1);endend5. 将温度随时间变化情况用曲线表示flnuiro604020.6 o6. 给出3000、9000、15000三个时刻的温度分布情况，对温度随时间变化规律做说明。T3000=100.0000 63.4362 34.2299 15.8021 7.464

8、1 7.4641 15.802134.2299 63.4362 100.0000T9000=100.0000 81.6930 65.6076 53.6839 47.3466 47.3466 53.683965.6076 81.6930 100.0000T15000=100.0000 89.9415 81.0962 74.5310 71.0378 71.0378 74.5310 81.0962 89.9415 100.0000根据数据分析，在同一个 x点上温度随时间的增加而增加，但增幅变小。x-T图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。7. 用计算数据说明，并结合差分方程余项，

9、空间、时间间隔对求解精度影响。数据量较大，且原理相同，我取一个向量演示一下。保持空间间隔不变，修改时间间隔，时间间隔加大，得到的误差加大。保持时间间隔不变，修改空间间隔，空间间隔加大，得到的误差加大。修改空间间隔的误差在增量比修改时间间隔的大。从方差余项上来看,Af d2u2 dt2(没有公式编辑器。只能从ppt里粘贴了)这个余项里的 t , x都在分母上，所以与误差成正比，且 x的次数 应该是比 t高，故影响较大。8. 用计算数据说明，稳定性要求对求解精度的影响。修改稳定性，即修改 x和t分的份数(ngrid ),之后看误差。稳定性越高，解的精度越高。即在满足稳定性要求 (a* t/( xA

10、2)<0.5 )时，a* t/( xA2)越接近0,误差越小。从概念上理解，稳定性越好，对引入时间层误差的抑制能力越强。所以误差越小。二、调用MATLAB函数完成上述计算1.编写M文件求解上述方程，并用适当的文字对程序做出说明。function pdepediaoyongm=0;x=linspace(0,1,11);%x 的网格t=linspace(0,20000,101);%t 的网格sol = pdepe(m,pdefun,icfun,bcfun,x,t);% 调用函数T=sol(:,:,1);% 解figure;% 画图surf(x,t,T)xlabel('x')y

11、label('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,:)%输出三个时间下的温度分布%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面function c,f,s=pdefun(x,t,T,DuDx)%PDE 方程函数c=100000;f=DuDx;s=0;function u0=ic

12、fun(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.将温度随时间变化情况用曲线表示。12010080604020024X 103.给出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.0000T9000=100.0000 83.4839 68.6032 56.8191 49.2705 46.6732 49.2705 56.8191 68.6032 83.4839 100.0000T15000=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图形仍为抛物 线，但随着时间的增加，极