BTCS.rar_btcs_有限差分_热传导_非定常_非定常热传导方程

dx=1/100; J1=99; x=(1:J1)*dx; xigema=0.1; dt=xigema*(dx^2); N=0.1/dt; t=(0:N-1)*dt; % % % % --------- ubc_l=0*ones(N,1); ubc_r=2*ones(N,1); for i=2:J1 if x(i)>0 & x(i)<0.3 f(i)=0; else if x(i)>=0.3 & x(i)<=0.6 f(i)=1; else f(i)=1+2.5*(x(i)-0.6); end end end u_initial=f(1,1:J1); % % % % --------- u=zeros(N,J1); u(1,:)=u_initial; % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % a:the main diagnal vector of coefficient matrix A;Au=f % b:the upper diagnal vector of A % c:the lower diagnal vector of A % f:Au=f a=ones(J1,1)*(2*xigema+1); b=-xigema*ones(J1-1,1); c=-xigema*ones(J1-1,1); f=zeros(J1,1); % % % % % % % % % for n1=2:N f=u(n1-1,:); f(1)=f(1)+xigema*ubc_l(n1); f(J1)=f(J1)+xigema*ubc_r(n1); u(n1,:)=qiujie(a,b,c,f);%追赶法求解三对角方程组 end [xx,tt]=meshgrid(x,t(1:5:end)); subplot(2,1,1) contourf(xx,tt,u(1:5:end,:)) colorbar title(('有限差分法(BTCS)')) k3=find(t==max(t)); T0=0; for k2=1:100 temp3=(2/k2/pi*cos(3*k2*pi/10)-5/k2/k2/pi/pi*sin(3*k2*pi/5))*exp(-k2^2*pi^2*k3)*sin(k2*pi*x); T0=T0+temp3; end T=T0+2*x; subplot(2,1,2) plot(x',T,'ro',x',u(k3,:),'*') mean_error=mean(abs(T-u(k3,:))); %平均误差 title(['t=',num2str(t(k3)),' mean_-error=',num2str(mean_error)]) %在标题中显示平均误差