计算热物理

目录
题目 (2)
1.方程和定解条件的无量纲化 (4)
2.控制方程的离散 (4)
3.稳定性分析 (5)
4.温度场的计算求解 (6)
5.温度场的等值线图 (9)
6.分析比较 (12)
7.个人心得体会 (14)
附录 (16)
题目：定义在正方形区域（如图所示）H*H=3m*3m 的二维非稳态导热问题的控制方程及其定解条件为
其中ρ，c，λ分别为导热体的密度、比热和导热系数，且ρ=7820kg/m3，c=460J/(kg•K)，λ=15W/(m•K)。

将定义域各方向均匀分成10个区块，取点中心网格分割方式，各边界共11个节点; 或者取块中心网格分割，各边界共10个节点，（两种网格，至少做一种），按以下要求，数值求解该问题。

（1）取（X，Y）=（x，y）
H ，θ=T−T0
T1−T0
，τ=t
ρcH2/λ
，将方程和定解条
件无量纲化，并在无量纲化条件下求解。

（2）采用以下方法对控制方程离散：
（a）用控制容积法离散成显式、全隐式格式；
（b）用ADI格式离散。

（3）用Fourier分析方法对以上3种格式的稳定性进行分析。

（4）嵌入初、边条件，分别对以上3种格式构成的定解问题编制计算程序，作数值求解，直到达到稳态温度分布——能直接求解显式格式和ADI格式构成的代数方程）直接求解，不能直接求解的（隐式格式构成的代数方程）采用以下方法迭代求解：
（a）Jacobi简单点迭代；
（b）Guess-Seidel点迭代；
（c）带松弛的Guess-Seidel点迭代；
（d）带松弛的Guess-Seidel线迭代；
（e）基于Guess-Seidel的交替方向线迭代。

并对收敛结果做出比较。

（5）选择全隐式格式在5个不同时刻（包括1个到达稳态的时刻）的计算结果，将无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。

从5个不同时刻的等值线图，分析时间发展中T的变化特性。

（6）根据数值计算情况，（精度高低，CPU时间消耗，编制程序难易等）对各种不同的离散方法及代数解法进行分析比较，并做出结论。

（7）综合上述内容和作业情况，写出作业报告，并对此类作业提出
建议和要求。

解：
（1）方程和定解条件的无量纲化：
将（X ，Y ）=（x ，y ）
H
，θ=
T−T 0T 1−T 0
，τ=
t
ρcH 2/λ
代入原方程组，无
量纲化，得到：
（2） 控制方程的离散： A ． 控制容积法
假定非稳态项中温度随空间为阶梯分布，扩散项中温度随时间线性分布、随空间阶梯分布 显式格式:
(θP −
θP 0)∆X∆Y
=(θE 0−θP
(
δX )e
−
θP 0−θW
(δX )w
)∆Y∆τ+(θN 0−θP 0
(δY )n
−
θP 0−θS
(δY )s
)∆X∆τ
写成：a P θP =a E θE 0+a W θW 0+a N θN 0+a S θS 0
+b 的形式，
则：a E =∆Y
(δX )e ，a W =∆Y
(
δX )w
，a N =∆X
(
δY )n
，a S =∆X
(
δY )s
a P =
∆X∆Y
∆τ
，a P 0=
∆X∆Y ∆τ
−a E −a W −a N −a S ，b =a P 0θP 0
全隐式格式：
(θP −θP 0)∆X∆Y =(θE
−θP (
δX )e
−
θP −θW (δX )w
)∆Y∆τ+(θN
−θ
P
(δY )
n
−θP −θS (δY )s
)∆X∆τ
写成：a P θP =a E θE +a W θW +a N θN +a S θS +b 的形式，
则：a E=∆Y
(δX)e ，a W=∆Y
(δX)w
，a N=∆X
(δY)n
，a S=∆X
(δY)s
a P=∆X∆Y
∆τ+a E+a W+a N+a S，a P0=∆X∆Y
∆τ
，b=a P0θP0
B．ADI格式：
θi,j∗−θi,j n ∆τ/2=
θi+1,j
∗−2θ
i,j
∗+θ
i−1,j
∗
∆X2
+
θi,j+1
n−2θ
i,j
n+θ
i,j−1
n
∆Y2
θi,j n+1−θi,j∗∆τ/2=
θi+1,j
∗−2θ
i,j
∗+θ
i−1,j
∗
∆X2
+
θi,j+1
n+1−2θ
i,j
n+1+θ
i,j−1
n+1
∆Y2
（3）稳定性分析:
A．显示格式：
a PθP=a EθE0+a WθW0+a NθN0+a SθS0+b
a P g=a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y+a P0(*)
本题中定义域各方向均匀分成10个区块，则
(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1
a E=a W=a N=a S=1，a P=∆X∆Y
∆τ
则(*)式可化为：g=2cos(k X∆X)+2cos(k Y∆Y)+(∆X∆Y/∆τ−4)
∆X∆Y/∆τ
假设k X=k Y，则g=1−8∆τ
∆X2sin2(k∆X
2
)
有|g|≤1恒成立，则∆τ
∆X2sin2(k∆X
2
)≤1
4
∆τ∆X2≤1
4
，∆τ≤1
4
∆X2。

即∆τ≤1
4
∆X2时，该显式格式收敛。

本题中，取∆X=0.1，故∆τ≤0.0025。

B．全隐式格式：
a PθP=a EθE+a WθW+a NθN+a SθS+b
a P g=g(a E e ik X∆X+a W e−ik X∆X+a N e ik Y∆Y+a S e−ik Y∆Y)+a P0(**)
本题中定义域各方向均匀分成10个区块，则
(δX)e=(δX)w=(δY)n=(δY)s=∆X=∆Y=0.1
a E=a W=a N=a S=1，a P0=∆X∆Y
∆τ
(**)式化为：g=∆X∆Y/∆τ
4−2cos(k X∆X)−2cos(k Y∆Y)+∆X∆Y/∆τ显然|g|≤1恒成立，故全隐式格式恒稳定。

C．ADI格式：
θi,j∗−θi,j n ∆τ/2=
θi+1,j
∗−2θ
i,j
∗+θ
i−1,j
∗
∆X2
+
θi,j+1
n−2θ
i,j
n+θ
i,j−1
n
∆Y2
θi,j n+1−θi,j∗∆τ/2=
θi+1,j
∗−2θ
i,j
∗+θ
i−1,j
∗
∆X2
+
θi,j+1
n+1−2θ
i,j
n+1+θ
i,j−1
n+1
∆Y2
空间导数离散式用微分算子符号写出，上面两个等式可化为：
(1−αX
2δX2)θi,j∗=(1+αY
2
δY2)θi,j n，αX=∆τ/∆X2
(1−αY
2δY2)θi,j n+1=(1+αX
2
δX2)θi,j∗，αY=∆τ/∆Y2
消去θi,j∗，得：
(1−αX
2
δX2)(1−
αY
2
δY2)θi,j n+1=(1+
αX
2
δX2)(1+
αY
2
δY2)θi,j n
故g=(1−2αX sin2k X∆X
2
)(1−2αY sin2k Y∆Y
2
)
(1+2αX sin2k X∆X
2
)(1+2αY sin2k Y∆Y
2
)
显然|g|≤1恒成立，故ADI格式恒稳定。

（4）温度场的计算求解
将定义域各方向均匀分成10 个区块,采用点中心网格分割方式，各边界共11 个节点。

边界条件的离散：
（1）Y=0：θi,1=1(i=1—11)；
（2）Y=1：θi,11=0(i=1—11)；
（3）X=0：j=2—10
显式：a 1,j θ1,j =a 2,j θ2,j +a 1,j+1θ1,j+1+a 1,j−1θ1,j−1+b 其中： a 1,j =
∆X∆Y 2∆τ
，a 2,j =∆Y
(
δX )e
，a 1,j+1=
∆X 2(δY )n
，a 1,j−1=
∆X 2(δY )s
b =θ1,j 0
(a 1,j −a 2,j −a 1,j+1−a 1,j−1)+q w ∆Y
隐式：a 1,j θ1,j =a 2,j θ2,j +a 1,j+1θ1,j+1+a 1,j−1θ1,j−1+b 其中： a 2,j =∆Y (
δX )e ，a 1,j+1=
∆X 2(δY )n
，a 1,j−1=
∆X 2(δY )s
a 1,j =
∆X∆Y
2∆τ
+ a 2,j +a 1,j+1+a 1,j−1，b =
∆X∆Y 2∆τ
θ1,j 0
+q w ∆Y
（4）X=1：j=2—10
显式：a 11,j θ11,j =a 10,j θ10,j +a 11,j+1θ11,j+1+a 11,j−1θ11,j−1+b 其中：a 11,j =
∆X∆Y 2∆τ
，a 10,j =∆Y
(
δX )w
，a 11,j+1=
∆X 2(δY )n
，a 11,j−1=
∆X 2(δY )s
b =θ11,j 0
(a 11,j −a 10,j −a 11,j+1−a 11,j−1)+q e ∆Y
隐式：a 11,j θ11,j =a 10,j θ10,j +a 11,j+1θ11,j+1+a 11,j−1θ11,j−1+b 其中： a 10,j =∆Y
(
δX )w
，a 11,j+1=
∆X 2(δY )n
，a 11,j−1=
∆X 2(δY )s
a 11,j =
∆X∆Y 2∆τ
+ a 10,j +a 11,j+1+a 11,j−1，b =
∆X∆Y 2∆τ
θ11,j 0
+q e ∆Y
初始条件：
由“初始T(x,y,0)为x=0→x=H 的温度成线性分布”我们可以推知：
θi,j =1−0.1∗(i −1) i,j=1—11
据此编程，编程软件：Visual C++，程序代码见附录。

其中，左右边界的热流量qe 和qw 、时间步长、稳态精度、迭代精度、时间上限以及最大迭代次数均可在define 项目中直接调整。

几种迭代算法的比较：
1.Jacobi简单点迭代：
利用上一时刻的节点温度值计算下一时刻该节点的温度，程序较为简单。

收敛速度慢，迭代最大次数为9次，出现在时间刚开始的阶段。

下面是计算结果：
2.Gauss-Seidel点迭代：
在计算节点值时，采用邻近节点的最新可用值，而不局限于前一迭代完成后所得到的值。

此种方法由于及时引入新值，迭代收敛速度快于简单迭代。

在实际程序设计中，先计算边界，使边界的最新可用值被尽可能快的引入到内节点的计算中，对加快收敛速度起到很好的作用。

在以后的程序中，也都采用了这种方法。

该算法的最大迭代次数为7次，出现在时间值很小的起步阶段。

下面是计算结果：
3.松弛迭代
本次上机取松弛因子w在从0.1开始按0.1增量递增到1.9。

松弛迭代的收敛速度与松弛因子有关。

当松弛因子在1.2-1.4时，收敛速度明显加快，而当松弛因子取0.1或1.9附近的值时，有可能发散。

4.Gauss-Seidel的交替方向线迭代
此法扫描方向可以变化，且有多种组合。

本程序先处理边界条件，然后采用从左至右的列扫描，后再采用从上至下的行扫描，全场两次扫描组成一轮迭代，这种扫描是对一条线上隐式格式迭代求解。

相对于Gauss-Seidel迭代，收敛速度有所提高，最大迭代次数为4次。

下面是计算结果：
（5）温度场的等值线图
这里，结合Gauss-Seidel的交替方向线迭代的计算结果，无量纲计算量还原成有量纲量，分别画出温度T在平面空间的等值线图。

我们可以看到，达到稳态所需的无量纲时间为0.636，折算成实际时间为38小时。

等值线图如下所示：
T=0时（初始时刻）：
T=7.6h时：
T=15.2h时：
T=22.8h时：
T=30.4h时：
T=38.0h时（稳态时刻）：
从上面的图可以看出：
（1）刚开始，同水平线温度相同，同一竖直线上温度呈线性分布，符合初始条件的要求。

（2）在时间间隔相同的情况下，刚刚开始时温度变化比较迅速，随着时间增加，温度变化趋向缓慢，最后达到稳定。

（3）由于左右边界条件对称，理论上温度分布也是左右对称。

数值解得到的结果也验证了这一点。

最高温度为1.0885，
折算成实际温度为413.275℃，出现在图的左下角和右下
角。

（4）随着时间的推进，低温区逐渐缩小，高温区逐渐扩大。

与左右边界有恒定热流进入此区域有关。

综上可见，计算结果与实际情况相符。

（6）分析比较
这里先把显示格式和ADI格式的计算结果予以展示，以作为对前面隐式格式迭代求解计算结果的对比。

显示格式计算结果：
ADI格式计算结果：
A．精度
从几种算法的运算结果来看，精度上的差别不大；但从达到稳态所需时间来看，精度上还是有差别。

计算结果的精度主要取决于离散方程的精度，由离散方程截断误差的量级所决定。

本题中所使用的各种迭代算法的精度均为二阶精度。

时间方面，由于求解下一时刻温度场时，显示格式和ADI格式上一时刻的温度场并没有收敛，而隐式格式上一时刻的温度场是收敛的，两者的差异导致最终达到稳态所需的时间不同，相差0.07，约为33min。

B．CPU时间消耗
显示格式和ADI格式的求解不需要迭代，运算量小，所
以耗时短，分别为0.062秒和0.078秒。

相比之下，隐式格式的算法普遍耗时长。

耗时最多的是Gauss-Seidel迭代，这是因为相比于Jacobi迭代，前者对边界的处理采用了一种更容易理解的方式，但需要额外的空间，也增加了运算量，所以虽然迭代次数有所减少，但总的运行时间反而增加。

Gauss-Seidel 线迭代的运行时间明显比Gauss-Seidel迭代减少很多，这主要归功于迭代次数的进一步减少。

C．编程难易程度
显示格式的程序设计起来简单，不需要考虑迭代等一系列问题。

ADI格式也是直接求解，程序较为简单。

隐式格式的求解，相对前两者就复杂了一些，各种算法间也有较大差异。

Jacobi迭代和Gauss-Seidel迭代算法都较简单，编程难度几乎同等，Gauss-Seidel线迭代的算法相对难一些，编程难度也大了一些。

而带松弛因子的迭代，需要在时间步长和迭代次数之上再增加一个变量，程序设计的难度大大增加，其中带松弛的Gauss-Seidel线迭代的程序部分是本人认为最难设计的。

结论：
显示和ADI格式，虽然可以直接求解，但精度不高，受时间步长的影响也很大。

带松弛的迭代需要考虑松弛因子的选取，若选取不当将会带来巨大的运算量，而且程序设计难度大。

隐式格式不受时间步长的影响，Gauss-Seidel线迭代的迭代次数少，CPU运行时间短，程序设计难度可以接受。

综
合考虑，本人认为隐式格式Gauss-Seidel线迭代是最为理想
的求解方法。

（7）个人心得体会
此次计算热物理大作业的完成，前后用了将近一周的时间，难度最大的部分在于程序源代码的编写。

在代码的编写工程中，我们需要考虑多种参量与变量，以及各种各样的算法。

其中，最难处理的部分在于边界条件，这不仅仅因为需要考虑的因素较内节点多，还因为需要编写专门的算法，因而增加了程序的设计难度，而且，稍有不慎就满盘皆输。

通过好几百行代码的编写，巩固了既有的程序设计基础，又学到了一些新的东西，如程序执行时间—CPU耗时的计算。

就计算热物理这门学科而言，通过这次大作业，将课堂所学加以实际应用，学以致用，不仅扎实了基础，还得到了一些感悟和收获。

最后，说一点个人的看法：二维数组的求解，最起码的时间复杂度为N的平方项，如果考虑迭代，增加迭代次数，则又增加了一阶复杂性，现在我们所处理的非稳态问题，本身就需要考虑时间的递增，至此，我们可以保守估计，程序的时间复杂度已经达到了N的4次方。

如果再加入可变的松弛因子，复杂度就达到5次方。

时间复杂度，其实也代表了程序设计的复杂性，一个5重循环就足以让人头疼不已，更何况还得考虑多种复杂变化的变量。

所以，我建议在大作业中舍去松弛因子这一项，使复杂度降阶。

如果要让我们认识到松弛因子的作用以及不同松弛因子对迭
代次数的影响，可以专门设计一个不包含时间项的问题。

附录：（源代码部分）
1.显式格式：
#include<iostream>
#include<iomanip>
#include<ctime>
#define dx 0.1 //x方向步长
#define dy 0.1 //y方向步长
#define qw 1.0
#define qe 1.0
#define n (int)(1/dx)+2 //n=12
#define tup 0.002 //时间步长，不大于0.0025
#define tup_limit 50 //时间上限
#define epsilon 0.00001
using namespace std;
//数组初始化
void init(double T[n][n]){
int i,j;
for(i=1;i<n;i++)
for(j=1;j<n;j++)
T[i][j]=1-0.1*(i-1);
}
//打印数组
void print(double T[n][n]){
int i,j;
for(i=n-1;i>0;i--){
for(j=1;j<n;j++)
cout<<setw(5)<<T[i][j]<<" ";
cout<<endl;
}
}
//绝对值函数
double re(double x){
if(x<0.0)x=0.0-x;
return x;
}
//向量范数
int fanshu(double X[],double Y[]){
int i;
double s=0.0;
for(i=1;i<n;i++)
if(s<re(X[i]-Y[i]))s=re(X[i]-Y[i]);
if(s<epsilon)return 0;
else return 1;
}
//计算温度值
void fun2(double T[n][n],double temp[n][n]){
int i,j;
double c1,c2;
c1=dx*dy/tup;
for(i=2;i<n-1;i++) //内部各点i,j=2-10
for(j=2;j<n-1;j++){
c2=temp[i][j-1]+temp[i][j+1]+temp[i-1][j]+temp[i+1][j];
c2+=temp[i][j]*(c1-4.0);
T[i][j]=c2/c1;
}
//左边界,TDMA方法求解
c1=c1/2.0;
temp[0][1]=0.0;
for(i=2;i<n-1;i++)temp[0][i]=0.5/(c1-0.5*temp[0][i-1]);
temp[1][2]=(temp[2][1]*(c1-2)+qw*dy+T[2][2]+0.5*T[1][1])/c1;
for(i=3;i<10;i++)temp[1][i]=(temp[i][1]*(c1-2)+qw*dy+T[i][2]+0.5*temp[1][i-1])/(c1-0.5*temp[0][i-1]);
T[10][1]=(temp[10][1]*(c1-2)+qw*dy+T[10][2]+0.5*T[11][1]+0.5*temp[1][9])/(c1-0.5*temp[0][9]);
for(i=9;i>1;i--)T[i][1]=temp[0][i]*T[i+1][1]+temp[1][i];
//右边界,TDMA方法求解
for(i=2;i<n-1;i++)temp[0][i]=0.5/(c1-0.5*temp[0][i-1]);
temp[1][2]=(temp[2][11]*(c1-2)+qe*dy+T[2][10]+0.5*T[1][11])/c1;
for(i=3;i<10;i++)temp[1][i]=(temp[i][11]*(c1-2)+qe*dy+T[i][10]+0.5*temp[1][i-1])/(c1-0.5*temp[0][i-1]);
T[10][11]=(temp[10][11]*(c1-2)+qe*dy+T[10][10]+0.5*T[11][11]+0.5*temp[1][9])/(c1-0.5*temp[0][9]);
for(i=9;i>1;i--)T[i][11]=temp[0][i]*T[i+1][11]+temp[1][i];
}
//时间向前走一步，并判断是否稳态
//与前一时刻的温度进行比较，若低于精度值则返回0，否则返回1
int fun1(double T[n][n]){
double temp[n][n];
int i,j;
for(i=1;i<n;i++)
for(j=1;j<n;j++)
temp[i][j]=T[i][j]; //保存前一时刻的值
fun2(T,temp); //计算该时刻的值
for(i=2;i<n-1;i++) //逐行比较温度值
if(fanshu(T[i],temp[i])!=0)return 1;
return 0;
}
//主函数
int main(){
clock_t start,end;
start = clock(); //程序运行时间
double T[n][n];
init(T); //数组初始化
//print(T); //打印数组
double time=0.0; //无量纲时间
while(time<=tup_limit){
time+=tup;
if(fun1(T)==0){
cout<<"无量纲时间"<<time<<"时达到稳态"<<endl;
cout<<"稳态时区域各节点的无量纲温度为:"<<endl;
cout<<setiosflags(ios::fixed)<<setprecision(4);
print(T);
end = clock();
cout<<"运行时间:"<<(double)(end - start) / CLOCKS_PER_SEC<<"S"<<endl;
return 0; //结束程序运行
}
}cout<<"在无量纲时间"<<tup_limit<<"内没有达到稳态."<<endl;
end = clock();
cout<<"运行时间"<<(double)(end - start) / CLOCKS_PER_SEC<<"秒"<<endl;
return 0;
}
2.ADI格式：
#include<iostream>
#include<iomanip>
#include<ctime>
#define dx 0.1 //x方向步长
#define dy 0.1 //y方向步长
#define qw 1.0
#define qe 1.0
#define n (int)(1/dx)+2 //n=12
#define tup 0.002 //时间步长，不大于0.0025 #define tup_limit 10 //时间上限
#define epsilon 0.00001
using namespace std;
//数组初始化
void init(double T[n][n]){
int i,j;
for(i=1;i<n;i++)
for(j=1;j<n;j++)
T[i][j]=1-0.1*(i-1);
}
//打印数组
void print(double T[n][n]){
int i,j;
for(i=n-1;i>0;i--){
for(j=1;j<n;j++)
cout<<setw(5)<<T[i][j]<<" ";
cout<<endl;
}
}
//绝对值函数
double re(double x){
if(x<0.0)x=0.0-x;
return x;
}
//向量范数
int fanshu(double X[],double Y[]){
int i;
double s=0.0;
for(i=1;i<n;i++)
if(s<re(X[i]-Y[i]))s=re(X[i]-Y[i]);
if(s<epsilon)return 0;
else return 1;
}
//计算温度值
void fun2(double T[n][n],double temp[n][n]){ int i,j;
double c1;
double p[9],q[9];
//左边界,TDMA方法求解
c1=dx*dy/tup;
p[0]=1.0/c1;
for(i=1;i<8;i++)p[i]=1.0/(c1-p[i-1]);
q[0]=(T[2][2]*2.0+2.0*qw*dy+T[1][1]+(c1-4)*temp[2][1])/c1;
for(i=1;i<8;i++)q[i]=(T[i+2][2]*2.0+2.0*qw*dy+(c1-4)*temp[i+2][1]+q[i-1])/(c1-p[i-1]);
T[10][1]=(T[10][2]*2.0+2.0*qw*dy+T[11][1]+(c1-4)*temp[10][1]+q[7])/(c1-p[7]);
for(i=9;i>1;i--)T[i][1]=p[i-2]*T[i+1][1]+q[i-2];
//右边界,TDMA方法求解
q[0]=(T[2][10]*2.0+2.0*qe*dy+T[1][11]+(c1-4)*temp[2][11])/c1;
for(i=1;i<8;i++)q[i]=(T[i+2][10]*2.0+2.0*qe*dy+(c1-4)*temp[i+2][11]+q[i-1])/(c1-p[i-1]);
T[10][11]=(T[10][10]*2.0+2.0*qe*dy+T[11][11]+(c1-4)*temp[10][11]+q[7])/(c1-p[7]);
for(i=9;i>1;i--)T[i][11]=p[i-2]*T[i+1][11]+q[i-2];
//内部各点ADI算法
c1=2.0*dx*dy/tup;
p[0]=1.0/(c1+2.0);
for(i=1;i<8;i++)p[i]=1.0/(c1+2.0-p[i-1]);
for(j=2;j<11;j++){ //列扫描
q[0]=(T[2][j]*(c1-2)+T[2][j+1]+T[2][j-1]+T[1][j])/(c1+2.0);
for(i=1;i<8;i++)q[i]=(T[i+2][j]*(c1-2)+T[i+2][j+1]+T[i+2][j-1]+q[i-1])/(c1+2.0-p[i-1]);
T[10][j]=(T[10][j]*(c1-2)+T[10][j+1]+T[10][j-1]+T[11][j]+q[7])/(c1+2.0-p[7]);
for(i=9;i>1;i--)T[i][j]=p[i-2]*T[i+1][j]+q[i-2];
}
//行扫描
for(i=2;i<11;i++){
q[0]=(T[i][2]*(c1-2)+T[i+1][2]+T[i-1][2]+T[i][1])/(c1+2.0);
for(j=1;j<8;j++)q[j]=(T[i][j+2]*(c1-2)+T[i+1][j+2]+T[i-1][j+2]+q[j-1])/(c1+2-p[j-1]);
T[i][10]=(T[i][10]*(c1-2)+T[i+1][10]+T[i-1][10]+T[i][11]+q[7])/(c1+2.0-p[7]);
for(j=9;j>1;j--)T[i][j]=p[j-2]*T[i][j+1]+q[j-2];
}
}
//时间向前走一步，并判断是否稳态
//与前一时刻的温度进行比较，若低于精度值则返回0，否则返回1
int fun1(double T[n][n]){
double temp[n][n];
int i,j;
for(i=1;i<n;i++)
for(j=1;j<n;j++)
temp[i][j]=T[i][j]; //保存前一时刻的值
fun2(T,temp); //计算该时刻的值
for(i=2;i<n-1;i++) //逐行比较温度值
if(fanshu(T[i],temp[i])!=0)return 1;。