传热学MATLAB温度分布大作业完整版
一维热传导MATLAB模拟
一维热传导MATLAB模拟昆明学院2015届毕业设计（论文）设计（论文）题目一维热传导问题的数值解法及其MATLAB模拟子课题题目无姓名伍有超学号 2所属系物理科学与技术系专业年级 2011级物理学2班指导教师王荣丽2015 年 5 月摘要本文介绍了利用分离变量法和有限差分法来求解一维传导问题的基本解，并对其物理意义进行了讨论。
从基本解可以看出，在温度平衡过程中，杠上各点均受初始状态的影响，而且基本解也满足归一化条件，表示在热传导过程中杆的总热量保持不变。
通过对一维杆热传导的分析，利用分离变量法和有限差分法对一维热传导进行求解，并用MATLAB 数学软件来对两种方法下的热传导过程进行模拟，通过对模拟所得三维图像进行取值分析，得出由分离变量法和有限差分法绘制的三维图基本相同，且均符合热传导过程中温度随时间、空间的变化规律，所以两种方法均可用来解决一维热传导过程中的温度变化问题。
关键词：一维热传导；分离变量法；有限差分法；数值计算；MATLAB 模拟AbstractIn this paper, the method of variable separation andfinite difference method are introduced to solve the problem of one-dimensional heat conduction problems, and the physical significance of numerical methods for heat conduction problems are discussed. From the basic solution, we can see the temperature on the bar are affected by the initial state during the process of temperature balance, and basic solution also satisfy the normalization condition which implied the invariance of the total heat in the bar during the heat conduction process. Through the analysis of the one-dimensional heat conduction, by taking use of variable separation method and finite difference method, we simulated the one-dimensional heat conduction problem by MATLAB. The three-dimensional images of the simulation results obtained by the method of separation of variables and finite difference method are similar to each other, and the temperature curve is in accordance with the law of temperature variation during heat conduction. Thus, we can go to the conclusion that both methods can be used to deal with the one-dimensional heat conduction problems.Keywords: One-dimensional heat conduction; method of variableseparation; finite difference method; numerical2method; MATLAB simulation目录第一章绪论11.1热传导的概念......................................................... .. (1)1.2热质的运动和传递......................................................... (1)第二章一维热传导问题的两种数值解法32.1一维热传导问题的初值问题32.2一维热传导问题的分离变量法42.3一维热传导问题的有限差分法63第三章一维有界杆热传导问题的MATLAB模拟9 3.1一维有界杆热传导问题93.2分离变量法的MATLAB模拟93.3有限差分法的MATLAB模拟12第四章总结与展望18参考文献19谢辞204第一章绪论1.1热传导的概念由于温度分布不均匀，热量从介质中温度高的地方流向温度低的地方称为热传导。
传热大作业 第4版4-23
东南大学能源与环境学院课程作业报告课程名称：传热学作业名称：传热学大作业——利用matlab程序解决热传导问题院（系）：能源与环境学院专业：热能与动力工程姓名：姜学号：完成时间：2012 年11 月8日评定成绩：审阅教师：目录一．题目及要求 (3)二．各节点离散化的代数方程..............................3&13 三．源程序......................................................5&16 四．不同初值时的温度分布情况...........................7&18 五．冷量损失的计算.......................................12&24 六．计算小结 (27)传热大作业——利用matlab 程序解决复杂热传导问题姓名:姜 学号: 班级:成绩：____________________一、题目及要求计算要求：一个长方形截面的冷空气通道的尺寸如附图所示。
假设在垂直于纸面的方向上冷空气及通道墙壁的温度变化很小，可以忽略。
试用数值方法计算下列两种情况下通道壁面中的温度分布及每米长度上通过壁面的冷量损失：（1） 内、外壁面分别维持在10℃及30℃；（2） 内、外壁面与流体发生对流传热，且有110f t C =︒、2120/()h W m K =⋅，230f t C =︒、224/()h W m K =⋅。
(取管道导热系数为0.025/()W m K λ=⋅)二、各节点的离散化的代数方程1、基本思想：将导热问题的温度场，用有限个离散点上的值的集合来代替，通过求解按一定方法建立起来的关于这些值的代数方程，来获得离散点上被求物理量的值。
2、基本步骤：(1)建立控制方程以及定解条件：对于(1)问有：2.2m3m 2m1.2m h 1、t f1h 1、t f2导热微分方程22220t t x y ∂∂+=∂∂定解条件为第一类边界条件对(2)问有： 导热微分方程22220t t x y ∂∂+=∂∂定解条件为第三类边界条件(2)区域离散化：如下图所示，用一系列与坐标轴平行的网格线把求解区域划分成许多子区域，以网格线的交点作为需要确定温度值的空间位置，称为节点。
方腔内自然对流MATLAB程序数值传热学
(2)
其中， 是热扩散系数。 边界条件： 由无滑移边界条件得，四个壁面上的速度均为零，即： un us uw ue 0 。 在热壁面上 T 100 ，在冷壁面上 T 0 ，在上下绝热壁面上处
T 0。 y
数值处理：
区域离散化如图 2 所示。 对于动量守恒方程，在不考虑压力的情况下先计算出一个临时速度
%========================================================================== uu(1:nx+1,1:ny+1)=0.5*(u(1:nx+1,1:ny+1)+u(1:nx+1,2:ny+2)); vv(1:nx+1,1:ny+1)=0.5*(v(1:nx+1,1:ny+1)+v(2:nx+2,1:ny+1));
%==========================================================================
un=0.0; us=0.0; vw=0.0; ve=0.0; Tw=100.0; Te=0.0; T(1,1:ny+1)=100; TT(1,1:ny+1)=100; %========================================================================== u(1:nx+1,1)=2*us-u(1:nx+1,2); u(1:nx+1,ny+2)=2*un-u(1:nx+1,ny+1); v(1,1:ny+1)=2*vw-v(2,1:ny+1); v(nx+2,1:ny+1)=2*ve-v(nx+1,1:ny+1);
(完整word版)matlab绘制温度场
通过在室内的某些位置布置适当的节点，采集回来室内的温湿度以及空气质量等实际参数。
首先对室内空间建模，用一个无限细化的三维矩阵来模拟出室内的温度分布情况，针对采集回来的数据，采用插值法和适当次数的拟合函数的拟合，得出三维矩阵的实际值的分布，最后结合matlab软件绘制出计算出的温度场的三维图像。
一．数据的采集与处理因为影响人的舒适感的温度层只是室内的某一高度范围内的温度，而温度传感器虽然是布置在一个平面内，但是采用插值法和拟合函数法是可以大致再现出影响人的舒适感的温度层的温度变化的。
同时，在构建出的三维模型中，用第三维表示传感器层面的温度。
在传感器层面，传感器分布矩阵如下：X=【7.5 36.5 65.5】（模型内单位为cm）Y=【5.5 32.5 59.5】Z=【z1 z2 z3；z4 z5 z6；z7 z8 z9；】（传感器采集到的实时参数）采用meshgrid（xi，yi，zi，…）产生网格矩阵；首先按照人的最小温度分辨值，将室内的分布矩阵按照同样的比例细化，均分，使取值点在坐标一定程度上也是接近于连续变化的，从而才能最大程度上使处理数据得来的分布值按最小分辨值连续变化！根据人体散热量计算公式：C=hc（tb-Ta）其中hc为对流交换系数；结合Gagge教授提出的TSENS热感觉指标可以计算出不同环境下人的对环境温度变化时人体温度感知分辨率，作为插值法的一个参考量，能使绘制出的温度场更加的符合人体的温度变化模式。
例如按照10cm的均差产生网格矩阵（实际上人对温度的分辨率是远远10cm大于这个值的，但是那样产生的网格矩阵也是异常庞大的，例如以0.5cm为例，那么就可以获得116*108=12528个元素，为方便说明现已10cm为例）：[xi yi]=meshgrid(7.5:10:65.5,5.5:10:59.5)xi =7.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.50007.5000 17.5000 27.5000 37.5000 47.5000 57.5000yi =5.5000 5.5000 5.5000 5.5000 5.5000 5.500015.5000 15.5000 15.5000 15.5000 15.5000 15.500025.5000 25.5000 25.5000 25.5000 25.5000 25.500035.5000 35.5000 35.5000 35.5000 35.5000 35.500045.5000 45.5000 45.5000 45.5000 45.5000 45.500055.5000 55.5000 55.5000 55.5000 55.5000 55.5000产生网格矩阵之后，就可以在测得的实时数据的基础上，通过相关的温度场的专业的估算函数，以及相关的数值处理函数来估计整个分布面（有最小的分辨率）上的温度了。
铜芯电缆温度分布MATLAB计算模型
1.题目如图1所示铜芯电缆，电流为5000A ，内径为10mm ，外包材料聚氯乙烯的厚度为2mm ，导热系数为0.15+0.00013{t}K)W/(m ⋅。
电缆左半边为绝热边界条件，右半边为第三类边界条件，空气温度为20℃，绝缘层表面与环境间的复合表面传热系数为10K)W/(m 2⋅。
铜的电阻率为()[]2010t-a R R t +=，m 1075.180⋅Ω⨯=-R ，C /004.0︒=a ，t 的单位为摄氏度。
试通过数值方法求解温度分布。
图 12.编程计算2.1 控制方程根据题意，本题为二维稳态导热问题，其控制方程为：011=+⎪⎭⎫ ⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂S T r r r T r r r θλθλ 边界条件：πππθ2~23,2~0007.0==r ：()W f B T T h q -=23~2007.0ππθ==r ：0=B q其中：C 20︒=f T 。
2.2 方程离散为建立通用方程，考虑非稳态项的控制方程为：011=+⎪⎭⎫⎝⎛∂∂∂∂+⎪⎭⎫ ⎝⎛∂∂∂∂=∂∂S T r r r T r r r t T cθλθλρ 采用全隐格式，在t ∆时间内，对控制容积积分，整理后可得：b T a T a T a T a T a S S N N W W E E P P ++++=其中：()ee e E r r a λδθ/∆=，()ww w W r r a λδθ/∆=，()n n n N r r a λδθ/∆=，()ss s S r r a λδθ/∆=，V S a a a a a a P P S N W E P ∆-++++=0()tVc a P ∆∆=ρ0，00P P C T a V S b +∆=，()θ∆∆+=∆r r r V s n 5.0采用通用表达式，各表达式如下表：表1 坐标及系数表达式2.3 边界条件处理对于北边界，采用附加源项法处理。
由于北边界(πππθ2~23,2~0007.0==r )为第三类边界条件，则最靠近边界的控制容积加入以下附加源项：()Bn f adC x h T V AS λδ//1,+∆=()Bn ad P x h V A S λδ//11,+∆-=其中：C 20︒=f T将附加源项加到相应控制容积后，再令相应的0N =a 。
matlab求解一维带内热源传热问题
matlab 求解一维带内热源传热问题解一维带有内部热源的传热问题通常涉及到热传导方程的求解。
热传导方程描述了温度场随时间和空间的变化。
一维热传导方程通常写作：22()T T Q x t xα∂∂=+∂∂ 其中：• T 是温度，• t 是时间，• x 是空间坐标，• α 是热扩散系数，• Q(x) 是热源。
解这个方程需要适当的边界条件和初始条件。
为了简化问题，我们可以考虑一个稳态（0T t∂=∂）情况。
以下是使用 MATLAB 求解一维带有内部热源的传热问题的简单示例代码：% 参数设置L = 1; % 区域长度alpha = 0.01; % 热扩散系数Q = @(x) 1; % 内部热源% 空间离散化N = 100; % 离散网格数x = linspace(0, L, N);% 热传导方程T = zeros(1, N);T(1) = 0; % 初始条件T(N) = 100; % 边界条件% 离散格式求解dx = x(2) - x(1);dt = 0.01;num_steps = 1000;for step = 1:num_stepsfor i = 2:N-1T(i) = T(i) + alpha * dt / dx^2 * (T(i+1) - 2*T(i) + T(i-1)) + Q(x(i)) * dt;endend% 结果可视化plot(x, T);xlabel('空间坐标');ylabel('温度');title('一维带内部热源传热问题');请注意，这是一个简化的例子，具体的问题可能需要更多的考虑，例如更精确的数值方法、不同的边界条件和初始条件、更复杂的热源分布等。
这个示例主要用于演示MATLAB 中解决这类问题的基本方法。
大作业题(2014版)——中科院传热学(雁栖湖校区)
（注意：如发现作业中互相拷贝及抄袭者，此次作业成绩全将记为 一矩形平板 0 x a , 0 y b ，内有均匀恒 定热源 g 0 ， 在 x 0 及 y 0 处绝热， 在x a及y b 处保持温度 T1 ， 初始时刻温度为 T0 ， 如右图 1 所示： 1、求 t 0 时，矩形区域内的温度分布 T x, y, t 的解 析表达式； 2、若 a 18m ， b 12m ， g 0 1W m3 ， T1 600 K ， T0 200 K ，热传导系数 热扩散系数 0.8 m 2 s 。 请根据 1 中所求温度分布用 MATLAB k 1.0W m K ， 软件绘出下列结果，并加以详细物理分析： (a) 300s 内，在同一图中画出点 (0,4) 、 (0,8) 、 6,0 、 (12,0) 、 (9,6) （单位：m） 温度随时间的变化 (b) 200s 内，画出点 (18,4) 、 (18,8) 、 6,12 、 (12,12) 、 (9,6) （单位：m）处，分 别沿 x、y 方向热流密度值随时间的变化 (c) 画出 t 50s、 75s、 100s、 125s、 150s 时刻区域内的等温线 (d) 300s 内，在同一图中画出点 9,0 ， 9, 6 （单位：m）在其他参数不变时 g 0 分 别等于 1W m 3 ， 2W m3 ， 3W m3 情况下的温度、热流密度 (e) 600s 内，在同一图中画出点(9,6)（单位：m）在其它参数不变情况时热导率 分别为 0.5W m K 、 1.0W m K 和 1.5W m K 时的温度、热流密度 (f) 600s 内，在同一图中画出点(9,6) （单位：m）在其它参数不变情况时热扩散 系数分别为 0.4 m 2 s 、 0.8 m 2 s 和 1.2 m 2 s 时的温度、热流密度 3、运用有限差分法计算 2 中(b)、(d)和(e)，并与解析解结果进行比较，且需将数 值解与解析解的相对误差减小到万分之一以下； 4、附上源！ ）
哈工程传热学数值计算大作业
传热学二维稳态导热问题的数值解法杨达文2011151419赵树明2011151427杨文晓2011151421吴鸿毅2011151416第一题：a=linspace(0,0.6,121);t1=[60+20*sin(pi*a/0.6)];t2=repmat(60,[80 121]);s=[t1;t2]; %构造矩阵for k=1:10000000 %理论最大迭代次数，想多大就设置多大S=s;for j=2:120for i=2:80S(i,j)=0.25*(S(i-1,j)+S(i+1,j)+S(i,j-1)+S(i,j+1));endendif norm(S-s)<0.0001break; %如果符合精度要求，提前结束迭代elses=S;endendS %输出数值解数值解数据量太大，这里就不打印出来，只画出温度分布。
画出温度分布：figure(1)xx=linspace(0,0.6,121);yy=linspace(0.4,0,81);[x,y]=meshgrid(xx,yy);surf(x,y,S)axis([0 0.6 0 0.4 60 80])grid onxlabel('L1')ylabel('L2')zlabel('t（温度）').60.66666777778L 1L 2t （温度）A0=[S(:,61)];for k=1:81B1(k)=A0(81-k+1);endB1 %x=L1/2时y方向的温度A1=[S(41,:)] %y=L2/2时x方向的温度x=0:0.005:0.6;y=0:0.005:0.4;A2=60+20*sin(pi*x/0.6)*((exp(pi*0.2/0.6)-exp(-pi*0.2/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6) )/2) %计算y=L2/2时x方向的解析温度B2=60+20*sin(pi*0.3/0.6)*((exp(pi*y/0.6)-exp(-pi*y/0.6))/2)/((exp(pi*0.4/0.6)-exp(-pi*0.4/0.6))/ 2) %计算x=L1/2时y方向的解析温度figure(2)subplot(2,2,1);plot(x,A1,'g-.',x,A2,'k:x'); %画出x=L1/2时y方向的温度场、画出x=L1/2时y方向的解析温度场曲线xlabel('L1');ylabel('t温度');title('y=L2/2');legend('数值解','解析解');subplot(2,2,2);plot(x,A1-A2); %画出具体温度场与解析温度场的差值曲线xlabel('L1');ylabel('差值');title('y=L2/2时，比较=数值解-解析解');subplot(2,2,3);plot(y,B1,'g-.',y,B2,'k:x'); %画出y=L2/2时x方向的温度场、画出y=L2/2时x方向的解析温度场曲线xlabel('L2');ylabel('t温度');title('x=L1/2');legend('数值解','解析解');subplot(2,2,4);plot(y,B1-B2); %画出具体温度场与解析温度场的差值曲线xlabel('L2');ylabel('差值');title('x=L1/2时，比较=数值解-解析解');y=L2/2时x方向的温度：60 60.1635347276130 60.3269574318083 60.4901561107239 60.653018915996160.8154342294146 60.9772907394204 61.1384775173935 61.298884093677961.4584005332920 61.6169175112734 61.7743263876045 61.930519281669662.0853891461909 62.2388298405943 62.3907362037523 62.541004126057762.6895306207746 62.8362138946214 62.9809534175351 63.123649991570263.2642058188844 63.4025245687647 63.5385114436490 63.672073244095163.8031184326565 63.9315571966177 64.0573015095482 64.180265191631864.3003639687311 64.4175155301449 64.5316395850212 64.642657917384664.7504944397430 64.8550752452343 64.9563286582797 65.054185283707565.1485780543131 65.2394422768254 65.3267156762441 65.410338438521565.4902532515567 65.5664053444751 65.6387425251668 65.707215216057165.7717764880854 65.8323820928694 65.8889904930310 65.941562890665265.9900632539310 66.0344583417471 66.0747177265744 66.110813815270166.1427218680003 66.1704200151959 66.1938892725421 66.213113553990066.2280796827826 66.2387774004857 66.2451993740203 66.247341200688866.2452014111934 66.2387814706441 66.2280857775556 66.213121660833566.1938993747528 66.1704320919304 66.1427358942990 66.110829762085766.0747355608048 66.0344780262737 65.9900847476605 65.941586148577365.8890154662295 65.8324087286383 65.7718047299493 65.707245003846265.6387737950858 65.5664380291767 65.4902872802189 65.410373736929465.3267521668755 65.2394798789402 65.1486166840471 65.054224854168964.9563690796505 64.8551164248743 64.7505362822981 64.642700324897664.5316824570463 64.4175587638655 64.3004074590802 64.180308831415964.0573451895733 63.9316008058186 63.8031618582281 63.672116371626463.5385541572596 63.4025667512431 63.2642473518283 63.123690755529062.9809932921539 62.8362527587866 62.6895683527611 62.541040603677462.3907713045038 62.2388634418130 62.0854211252013 61.930549515936761.7743547548873 61.6169438897778 61.4584248018242 61.298906131798361.1384972055701 60.9773079591820 60.8154488635041 60.653030848523060.4901652273162 60.3269636197632 60.1635378760476 60x=L1/2时y方向的温度：60 60.1308958471008 60.2618814819943 60.3930468323419 60.524481948785060.6562770664196 60.7885226663977 60.9213095376979 61.054728839108661.1888721614654 61.3238315901874 61.4596997681540 61.596569958966661.7345361106384 61.8736929197574 62.0141358961654 62.155961428198162.2992668485325 62.4441505006859 62.5907118062120 62.739051332642462.8892708622179 63.0414734614594 63.1957635516239 63.352246980097063.5110310927684 63.6722248074423 63.8359386883315 64.002285021688564.1713778926236 64.3433332631650 64.5182690516120 64.696305213238964.8775638224022 65.0621691561100 65.2502477791090 65.441928630549065.6373431122839 65.8366251788694 66.0399114293203 66.247341200688866.4590566635297 66.6752029193167 66.8959280998773 67.121383468913967.3517235256817 67.5871061108928 67.8276925149213 68.073647588380968.3251398551535 68.5823416279436 68.8454291264398 69.114582598162569.3899864420822 69.6718293350911 69.9603043614169 70.255609145064670.5579459853794 70.8675219958221 71.1845492460516 71.509244907413471.8418314019312 72.1825365549057 72.5315937512233 72.889242095483173.2557265760494 73.6312982331452 74.0162143310978 74.410738534857774.8151410909089 75.2296990126956 75.6546962706925 76.090423987246276.5371806363247 76.9952722483076 77.4650126199600 77.946723529732178.4407349585321 78.9473853161230 79.4670216732992 8066666666L 1t 温度y =L 2/2--1.--0.-3L 1差值y =L 2/2时，比较=数值解-解析解66778L 2t 温度x =L 1/200.050.10.150.20.250.30.350.4--1.--0.-3L 2差值x =L 1/2时，比较=数值解-解析解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
传热学大作业（第四章）姓名：张宝琪学号：03110608一、题目及要求1.各节点的离散化的代数方程2.源程序3.不同初值时的收敛快慢4.上下边界的热流量（λ=1W/(m℃））5.计算结果的等温线图6.计算小结题目：已知条件如下图所示：二、方程及程序(1)各温度节点的代数方程ta=(300+b+e)/4 ; tb=(200+a+c+f)/4; tc=(200+b+d+g)/4; td=(2*c+200+h)/4 te=(100+a+f+i)/4; tf=(b+e+g+j)/4; tg=(c+f+h+k)/4 ; th=(2*g+d+l)/4ti=(100+e+m+j)/4; tj=(f+i+k+n)/4; tk=(g+j+l+o)/4; tl=(2*k+h+q)/4tm=(2*i+300+n)/24; tn=(2*j+m+p+200)/24; to=(2*k+p+n+200)/24; tp=(l+o+100)/12 (2)源程序【G-S迭代程序】【方法一】函数文件为：function [y,n]=gauseidel(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);G=(D-L)\U;f=(D-L)\b;y=G*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=G*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0;0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0;0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0;0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0;0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]';[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6) xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.1523 84.1429 67.9096 63.3793 62.4214 20.1557 15.4521 14.8744 14.7746 【方法2】>> t=zeros(5,5);t(1,1)=100;t(1,2)=100;t(1,3)=100;t(1,4)=100;t(1,5)=100;t(2,1)=200;t(3,1)=200;t(4,1)=200;t(5,1)=200;for i=1:10t(2,2)=(300+t(3,2)+t(2,3))/4 ;t(3,2)=(200+t(2,2)+t(4,2)+t(3,3))/4;t(4,2)=(200+t(3,2)+t(5,2)+t(4,3))/4;t(5,2)=(2*t(4,2)+200+t(5,3))/4;t(2,3)=(100+t(2,2)+t(3,3)+t(2,4))/4;t(3,3)=(t(3,2)+t(2,3)+t(4,3)+t(3,4))/4; t(4,3)=(t(4,2)+t(3,3)+t(5,3)+t(4,4))/4; t(5,3)=(2*t(4,3)+t(5,2)+t(5,4))/4;t(2,4)=(100+t(2,3)+t(2,5)+t(3,4))/4;t(3,4)=(t(3,3)+t(2,4)+t(4,4)+t(3,5))/4;t(4,4)=(t(4,3)+t(4,5)+t(3,4)+t(5,4))/4;t(5,4)=(2*t(4,4)+t(5,3)+t(5,5))/4;t(2,5)=(2*t(2,4)+300+t(3,5))/24;t(3,5)=(2*t(3,4)+t(2,5)+t(4,5)+200)/24;t(4,5)=(2*t(4,4)+t(3,5)+t(5,5)+200)/24;t(5,5)=(t(5,4)+t(4,5)+100)/12;t'endcontour(t',50);ans =100.0000 200.0000 200.0000 200.0000 200.0000 100.0000 136.8905 146.9674 149.8587 150.7444 100.0000 102.3012 103.2880 103.8632 104.3496 100.0000 70.6264 61.9465 59.8018 59.6008 100.0000 19.0033 14.8903 14.5393 14.5117【Jacobi迭代程序】函数文件为：function [y,n]=jacobi(A,b,x0,eps)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);f=D\b;y=B*x0+f;n=1;while norm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end命令文件为：A=[4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0,0;-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0,0; 0,-1,4,-1,0,0,-1,0,0,0,0,0,0,0,0,0; 0,0,-2,4,0,0,0,-1,0,0,0,0,0,0,0,0;-1,0,0,0,4,-1,0,0,-1,0,0,0,0,0,0,0; 0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0,0; 0,0,-1,0,0,-1,4,-1,0,0,-1,0,0,0,0,0;0,0,0,-1,0,0,-2,4,0,0,0,-1,0,0,0,0;0,0,0,0,-1,0,-1,0,4,0,0,0,-1,0,0,0;0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0,0;0,0,0,0,0,0,-1,0,0,-1,4,-1,0,0,-1,0;0,0,0,0,0,0,0,-1,0,0,-2,4,0,0,0,-1;0,0,0,0,0,0,0,0,-2,0,0,0,24,-1,0,0;0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1,0;0,0,0,0,0,0,0,0,0,0,-2,0,0,-1,24,-1;0,0,0,0,0,0,0,0,0,0,0,-1,0,0,-1,12];b=[300,200,200,200,100,0,0,0,100,0,0,0,300,200,200,100]'; [x,n]=jacobi(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',1.0e-6); xx=1:1:4;yy=xx;[X,Y]=meshgrid(xx,yy);Z=reshape(x,4,4);Z=Z'contour(X,Y,Z,30)n =97Z =139.6088 150.3312 153.0517 153.5639108.1040 108.6641 108.3119 108.152384.1429 67.9096 63.3793 62.421420.1557 15.4521 14.8744 14.7746三、不同初值时的收敛快慢1、[方法1]在Gauss 迭代和Jacobi 迭代中，本程序应用的收敛条件均为norm(y-x0)>=eps ，即使前后所求误差达到e 的-6次方时，跳出循环得出结果。
将误差改为0.01时，只需迭代25次，如下[x,n]=gauseidel(A,b,[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]',0.01)运行结果为 将误差改为0.1时，需迭代20次，可见随着迭代次数增加，误差减小，变化速度也在减小。
[方法2]通过 i=1:10判断收敛，为迭代10次，若改为1:20，则迭代20次。
2、在同样的误差要求下，误差控制在e 的-6次方内，Gauss 迭代用了49次达到要求，而Jacobi 迭代用了97次，可见，在迭代中尽量采用最新值，可以大幅度的减少迭代次数，迭代过程收敛快一些。
在Gauss 中，初值为100，迭代46次达到精确度1.0e-6,初值为50时，迭代47次，初值为0时，迭代49次，初值为200时迭代50次，可见存在一个最佳初始值，是迭代最快。
这一点在jacobi 迭代中表现的尤为明显。
四、上下边界热流量：上边界t=200℃,∞t =10℃，所以，热流量Φ1=λ*[2*100-200x y∆∆+x ya∆∆t -200+x y∆∆bt -200+x y∆∆ct -200+2*t -200dx y∆∆]=1*（100/2+(200-139.6088)+(200-150.3312)+(200-153.0517)+(200-153.5639)/2) =230.2264W 下边界热流量Φ2=|λ*[x y∆∆m i t -t +x y∆∆o j t -t +x y∆∆p k t -t +2*t -t q l x y∆∆]-h*(2*10-100x y∆∆+x *t -t n ∆∆∞y+x *t -t o ∆∆∞y+x *t -t m ∆∆∞y+2*t -t p x y∆∆∞)|=|1*((84.1429-20.1557)+(67.9096-15.4521)+(63.3793-14.8744)+(62.4214-14.7746)/2)-10*(90/2+(20.1557-10)+(15.4521-10)+(14.8744 -10)+(14.7746-10)/2)| = |-489.925|W =489.25W五、温度等值线Gauss:Yacobi:六、计算小结导热问题进行有限差分数值计算的基本思想是把在时间、空间上连续的温度场用有限个离散点温度的集合来代替，即有限点代替无限点，通过求解根据傅里叶定律和能量守恒两大法则建立关于控制面内这些节点温度值的代数方程，获得各个离散点上的温度值。