1维非稳态导热隐式求解的matlab代码

clc;

clear;

tic

format longG

rho=100;

cp=1000;

k=10;

T0=3;

TW=3;

TE=5;

s=10;

t=60000;

x=3;

dt=100;

dx=0.1;

sizet=t/dt;

sizex=x/dx;

Tmatrix=zeros(sizet,sizex);

for i=1

   Tmatrix(i,:)=T0;

end

%我想用TDMA解方程组，所以需要一堆东西去记录ai1

%ap0始终不变，可以作为源项之一纳入b向量

%ap1、ae1和aw1不随时间变化但随坐标变化，所以应该以向量形式给出：

%ap1=zeros(1,)

%隐格式算一算

for i=2:sizet

   %设置一个sizex*sizex的系数矩阵A，用于存放ai1

   A=zeros(sizex,sizex);

   %结果向量b

   b=zeros(sizex,1);

   for j=1:1

       aw1=2*k/dx;%均分网格，dxe=dxw=dx;

       ae1=k/dx;

       ap0=rho*cp*dx/dt;

       ap1=ap0+ae1+aw1;

       A(j,j)=ap1;

       A(j,j+1)=-ae1;

       b(j)=s*dx+ap0*Tmatrix(i-1,j)+aw1*TW;

   end

   for j=2:sizex-1

       ap0=rho*cp*dx/dt;

       aw1=k/dx;

       ae1=k/dx;

       ap1=ap0+aw1+ae1;

       A(j,j-1)=-aw1;

       A(j,j)=ap1;

       A(j,j+1)=-ae1;

       b(j)=s*dx+ap0*Tmatrix(i-1,j);

   end

   for j=sizex:sizex

       ap0=rho*cp*dx/dt;

       aw1=k/dx;

       ae1=2*k/dx;

       ap1=ap0+aw1+ae1;

       %ap1*Tmatrix(i,j)-aw1*Tmatrix(i,j-1)=ap0*Tmatrix(i-1,j)+s*dx+ae1*TE;

       A(j,j-1)=-aw1;

       A(j,j)=ap1;

       b(j)=s*dx+ap0*Tmatrix(i-1,j)+ae1*TE;

   end

   A;

   T1=A\b;

   Tmatrix(i,:)=T1';

end

   Ans=Tmatrix(sizet,:)'

   toc