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

第五次作业（前三题写在作业纸上）

一、用有限差分方法求解一维非定常热传导方程，初始条件和边界条件见说明.pdf文件，热扩散系数α=const，

?T?2T??2 ?t?x1. 用Tylaor展开法推导出FTCS格式的差分方程

2. 讨论该方程的相容性和稳定性，并说明稳定性要求对求解差分方程的影响。 3. 说明该方程的类型和定解条件，如何在程序中实现这些定解条件。

4. 编写M文件求解上述方程，并用适当的文字对程序做出说明。（部分由网络搜索得

到，添加，修改后得到。） function rechuandaopde

%以下所用数据，除了t的范围我根据题目要求取到了20000，其余均从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%输出温度矩阵

dt=tspan(2)/ngrid(1);%t步长 h3000=3000/dt;

h9000=9000/dt;

h15000=15000/dt;000,9000,15000下，温度分别在T矩阵的哪些行 T3000=T(h3000,:) T9000=T(h9000,:)

T15000=T(h15000,:)%输出三个时间下的温度分布

%不再对三个时间下的温度-长度曲线画图，其图像就是三维图的截面

%稳定性讨论,傅里叶级数法 dx=xspan(2)/ngrid(2);%x步长 sta=4*a*dt/(dx^2)*(sin(pi/2))^2; if sta>0,sta<2

fprintf('\\n%s\\n','有稳定性') else

fprintf('\\n%s\\n','没有稳定性') 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,Te,jmax) for j=1:jmax

rms=((T(j)-Te(j))^2/jmax)^(1/2) end

function[Ue,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)=100-(400/(2*g-1)/pi)*sin((2*g-1)*pi*xe(j))*exp(-a*(2*g-1)^2*pi^2*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=zeros(ngrid); U(:,1)=g1(t);%边界条件 U(:,m)=g2(t);

U(1,:)=f(x);%初值条件 %差分计算 for j=2:n for i=2:m-1

U(j,i)=(1-2*a*k/h^2)*U(j-1,i)+a*k/h^2*U(j-1,i-1)+a*k/h^2*U(j-1,i+1); end end

5. 将温度随时间变化情况用曲线表示

1008060T4020021.5x 104110.5t000.20.4x0.60.8

6. 给出3000、9000、15000三个时刻的温度分布情况，对温度随时间变化规律做说明。 T3000=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图形仍为抛物线，但随着时间的增加，极值变小，图像变得平缓。

7. 用计算数据说明，并结合差分方程余项，空间、时间间隔对求解精度影响。 数据量较大，且原理相同，我取一个向量演示一下。

保持空间间隔不变，修改时间间隔，时间间隔加大，得到的误差加大。 保持时间间隔不变，修改空间间隔，空间间隔加大，得到的误差加大。 修改空间间隔的误差在增量比修改时间间隔的大。

从方差余项上来看，（没有公式编辑器。。。。只能从

ppt里粘贴了）这个余项里的△t，△x都在分母上，所以与误差成正比，且△x的次数应该是比△t高，故影响较大。

8. 用计算数据说明，稳定性要求对求解精度的影响。 修改稳定性，即修改x和t分的份数（ngrid），之后看误差。

稳定性越高，解的精度越高。即在满足稳定性要求（a*△t/(△x^2)<0.5）时，a*△t/(△x^2)越接近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(:,:,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;000,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=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. 将温度随时间变化情况用曲线表示。

12010080604020021.5x 104110.5000.20.40.60.8

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.0000

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函数的稳定性要用什么检验。傅里叶级数检验不适用。

本文来源：https://www.bwwdw.com/article/2v4r.html