一维稳态导热问题数值模拟

一维稳态导热问题数值模拟问题描述：设有一导热方程，边界条件为 编写一段程序对此问题进行数值模拟。解析： 1、用控制容积有限差分方法做出内部节点和边界节点的离散化方程：首先进行离散化，先确定节点，再确定控制容积。将0-1划分为N段，共N+1个节点，N个控制容积，其中。可以得到如下：对原方程建立差分方程，内部节点有：则转换为下式，： i=2，.,N上式即为内部节点的离散化方程。对于外部节点可有：综上可以得到内部节点和外部节点的离散化方程为：即为上式不满足系数为负数，则可改用如下离散方程：内部节点： 边界节点 边界节点 组成代数方程组： 写成矩阵方程组：2、写出代数方程组的迭代求解程序:用Matlab编写如下求解程序；（1）高斯赛德尔迭代法（调用程序gauseidel文件）function x,n=gauseidel(A,b,x0,eps,M)%高斯迭代格式%线性方程组的系数：A%线性方程组中常数向量：b%迭代初始向量：x0%解的精度控制：eps%迭代步数控制：M%线性方程组的解：x%求出所需精度的解实际迭代步数：nif nargin=3 eps=0.; M=10000;elseif nargin=4 M=10000;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);G=(D-L)U;f=(D-L)b;x=x0;n=0;tol=1;while tol=eps x=G*x0+f; n=n+1; tol=norm(x-x0); x0=x; if (n=M) disp (Warning：迭代次数过多，可能不收敛.) return; endend（2）主程序（demo文件）如下：N=input(请输入N值n)Tp=input(请输入Tp值n)x1=zeros(N,1)A0=zeros(N);A0(1,1)=N+1/(2*N);A0(1,2)=-N;A0(N,N-1)=-N;A0(N,N)=2*N+1/N;for i=2:N-1 A0(i,i-1)=-N; A0(i,i)=2*N+1/N; A0(i,i+1)=-N;endb0=zeros(N,1);b0(1,1)=(1/N)*Tp;b0(N,1)=(2/N)*Tp+N;for i=2:N-1 b0(i,1)=(2/N)*Tp;endA=A0;b=b0;x0=x1;x,n=gauseidel(A,b,x0)x=x;1t=(0:1/N:1)title(一维稳态导热问题空间温度分布图)xlabel(空间分布X)ylabel(温度分布T)hold onplot(t,x) 1.3517 1.3485 1.3387 1.3224 1.2992 1.2691 1.2316 1.1865 1.1332 1.0712方程组的解3、结果分析，以上程序计算当取=1。 当N=10时，迭代次数n=343 1.3519 1.3510 1.3486 1.3446 1.3389 1.3315 1.3225 1.3118 1.2993 1.2851 1.2692 1.2514 1.2317 1.2101 1.1865 1.1609 1.1332 1.1034 1.0713 1.0368方程组的解当N=20时，迭代次数n=1243 1.3519 1.3515 1.3504 1.3486 1.3461 1.3428 1.3389 1.3341 1.3287 1.3225 1.3155 1.3078 1.2993 1.2901 1.2800 1.2692 1.2575 1.2450 1.2317 1.2175 1.2025 1.1865 1.1697 1.1519 1.1332 1.1136 1.0929 1.0713 1.0486 1.0248方程组的解 当N=30时，迭代次数n=2634 以上图形是不同步长