[数值分析]上机实验报告
数值分析上机实验报告
《数值分析》上机实验报告
1. 用Newton 法求方程 X 7-X 4+14=0
在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据：
设函数在有限区间[a，b]上二阶导数存在，且满足条件
1. f (x ) f (b ) >0
2. f .. (x ) 在区间[a , b ]上不变号3. f . (x ) ≠0;
|f (c ) |
≤|f . (x ) |,其中c 是a , b 中使m ir (|f . (a ), f . (b ) |)达到的一个b -a
则对任意初始近似值x 0∈[a , b ],由Newton 迭代过程4.
x k +1=ϕ(x k ) =x k -
f (x k )
, k =0, 1, 2, 3 f ' (x k )
{x k }平方收敛于方程f (x ) =0在区间[a , b ]上的惟一解α所生的迭代序列
令
f (x ) =x 7-28x 4+14, f (0. 1) >0, f (1. 9) 0
5
2
2
3
故以1.9为起点
f (x k ) ⎧
⎪x k +1=x k -
f '(x k ) ⎨
⎪x 0=1. 9⎩
如此一次一次的迭代，逼近x 的真实根。当前后两个的差
1.2 C 语言程序原代码：
#include #include main()
{double x2,f,f1;
double x1=1.9; //取初值为 1.9 do
{x2=x1;
f=pow(x2,7)-28*pow(x2,4)+14; f1=7*pow(x2,6)-4*28*pow(x2,3); x1=x2-f/f1;}
while(fabs(x1-x2)>=0.00001||x1
1.3 运行结果：
1.4 MATLAB 上机程序
function y=Newton(f,df,x0,eps,M) d=0;
for k=1:M
if feval(df,x0)==0 d=2;break else
x1=x0-feval(f,x0)/feval(df,x0); end
e=abs(x1-x0); x0=x1;
if e
if d==1 y=x1; elseif d==0
y=' 迭代M 次失败' ; else
y= ' 奇异' end
function y=df(x) y=7*x^6-28*4*x^3; End
function y=f(x) y=x^7-28*x^4+14; End
>> x0=1.9;
>> eps=0.00001; >> M=100;
>> x=Newton('f','df',x0,eps,M); >> vpa(x,7)
1.5 问题讨论：
1. 使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。此程序的不足之处是，所要求解的方程必须满足上述定理的四个条件，但是第二和第四个条件在计算机上比较难以实现。
2.Newton 迭代法是一个二阶收敛迭代式，他的几何意义Xi+1是Xi 的切线与x 轴的交点, 故也称为切线法。它是平方收敛的，但它是局部收敛的，即要求初始值与方程的根充分接近，所以在计算过程中需要先确定初始值。
3. 本题在理论依据部分，讨论了区间(0.1,1.9)两端点是否能作为Newton 迭代的初值，结果发现0.1不满足条件，而1.9满足，能作为初值。另外，该程序简单，只有一个循环，且为顺序结构，故采用do-while 循环。当然也可以选择for 和while 循环。
2. 已知函数值如下表：
试用三次样条插值求f(4.563)及f ’(4.563)的近似值。 2.1 理论依据
S (x ) =M j -1
(x j -x ) 36h j -1
+M j
(x -x j -1) 3
6h j -1
h 2x j -x h 2x -x j -1j -1j -1
+(y j -1-M j -1)() +(y j -M j )()
6h j -16h j -1
这里h j -1=x j -x j -1 ，所以只要求出M j ，就能得出插值函数S （x ）。
⎡21
⎢μ
⎢12λ1⎢μ22
求M j 的方法为：⎢
⎢ ⎢⎢⎢⎣
⎤⎡M 0⎤⎡d 0⎤
⎥⎢M ⎥⎢d ⎥⎥⎢1⎥⎢1⎥⎥⎢⎥⎢⎥⎥⎢⎥=⎢⎥
⎥⎢⎥⎢⎥
⎥⎢⎥2λN -1⎥⎢⎥⎢⎥⎢⎥M 12⎥⎢⎦⎣N ⎥⎦⎢⎣d N ⎥⎦
λ2
μN -1
6y 1-y 0⎧
') d =-y 0
⎪0h (h
00
⎪⎪y j +1-y j y j -y j -16d =(-) (j =1, 2, , N -1) j
这里⎪ h +h h h j -1j j j -1⎪
⎨
⎪d =6[y '-1(y -y )]
N -1
⎪N h N -1N h N -1N ⎪
h j
⎪μ=h j -1
λ=1-μ=j j ⎪j h +h h j -1+h j j -1j ⎩
最终归结为求解一个三对角阵的解。
用追赶法解三对角阵的方法如下：
⎤⎡b 1c 1⎤⎡1⎤⎡β1γ1
⎥⎢a b ⎥⎢l 1⎥⎢c γβ222122⎥⎢⎥⎢⎥⎢
⎥=LU ⎥=⎢⎥⎢ A =⎢l 21
⎥⎢⎥⎢⎥⎢
b c a βγn -1n -1⎥n -1n -1n -1⎥⎢⎢⎥⎢
⎢l n 1⎥a n b n ⎥βn ⎥⎣⎦⎢⎣⎦⎢⎣⎦
⎡δ1⎤
⎧L δ=d ⎥, 则由L δ=d 得 LUx =d , 即⎨, 若记δ=⎢ ⎢⎥⎩Ux =δ⎢⎣δn ⎥⎦
⎡1⎤⎡δ1⎤⎡d 1⎤
⎢l ⎥⎢ ⎥⎢ ⎥12⎢⎥⎢⎥=⎢⎥ ， ⎢ ⎥⎢ ⎥⎢ ⎥⎢⎥⎢⎥⎢⎥
l 1n ⎣⎦⎣δn ⎦⎣d n ⎦⎡β1
⎢⎢⎢⎢⎣
γ1
βn -1
⎤⎡x 1⎤⎡δ1⎤⎥⎢ ⎥⎢ ⎥⎥⎢⎥=⎢⎥ γn -1⎥⎢ ⎥⎢ ⎥
⎥⎢⎥⎢⎥βn ⎦⎣x n ⎦⎣δn ⎦
综上可得求解方程Ax=d的算法：
αi +1⎧
β=b , δ=d , l =, βi +1=b i +1-l i +1c i , 11i +1⎪11
βi
⎪⎪
i =1, 2,3, , n -1 ⎨δi +1=d i +1-l i +1δi
⎪δδ-c x
⎪x n =n , x i =i i i +1, i =n -1, , 2,1
βn βi ⎪⎩
2.2 C语言程序代码：
#include
#include
void main() {int i,j,m,n,k,p;
double q10,p10,s4,g4,x0,x1,g0=1,g9=0.1;; double s[10][10];
double a[10],b[10],c[10],d[10],e[10],x[10],h[9],u[9],r[9];
double f[10]={0,0.69314718,1.0986123,1.3862944,1.6094378, 1.7917595,1.9459101,2.079445,2.1972246,2.3025851}; printf("请依次输入xi:\n"); for(i=0;i
scanf("%lf",&e[i]); //求h 矩阵 for(n=0;n
d[0]=6*((f[1]-f[0])/h[0]-g0)/h[0];
d[9]=6*(g9-(f[9]-f[8])/h[8])/h[8]; for(j=0;j
d[j+1]=6*((f[j+2]-f[j+1])/h[j+1]-(f[j+1]-f[j])/h[j])/(h[j]+h[j+1]); for(m=1;m
u[m]=h[m-1]/(h[m-1]+h[m]); for(k=1;k
for(i=0;i
if(i==p)s[i][p]=2;} s[0][1]=1; s[9][8]=1; for(i=1;i
printf("三对角矩阵为:\n"); for(i=0;i
for(p=0;p
{printf("\n");} }
printf("根据追赶法解三对角矩阵得：\n"); a[0]=s[0][0]; b[0]=d[0]; for(i=1;i
{c[i]=s[i][i-1]/a[i-1]; //求d 矩阵 a[i]=s[i][i]-s[i-1][i]*c[i]; b[i]=d[i]-c[i]*b[i-1]; if(i==8) {p10=b[i]; q10=a[i];}} x[9]=p10/q10;
printf("M[10]=%lf\n",x[9]); for(i=9;i>=1;i--)
{x[i-1]=(b[i-1]-s[i-1][i]*x[i])/a[i-1]; printf("M[%d]=%lf\n",i,x[i-1]);}
printf("可得s(x)在区间[4,5]上的表达式;\n"); printf("将x=4.563代入得：\n"); x0=5-4.563; x1=4.563-4;
s4=x[3]*pow(x0,3)/6+x[4]*pow(x1,3)/6+(f[3]-x[3]/6)*(5-4.563)+(f[4]-x[4]/6)*(4.563-4);
g4=-x[3]*pow(x0,2)/2+x[4]*pow(x1,2)/2-(f[3]-x[3]/6)+(f[4]-x[4]/6);
printf("计算结果：f(4.563)的函数值是：%lf\nf(4.563)的导数值是：%lf\n",s4,g4);}
2.3 运行结果：
2.4 问题讨论
1. 三次样条插值效果比Lagrange 插值好，没有Runge 现象，光滑性较好。 2. 本题的对任意划分的三弯矩插值法可以解决非等距节点的一般性问题。 3. 编程过程中由于定义的数组比较多，需要仔细弄清楚各数组所代表的参数，要注意各下标代表的含义，特别是在用追赶法计算的过程中。
3. 用Romberg 算法求⎰13x x 1. 4(5x +7) sin x 2dx (允许误差ε=0. 00001) . 3.1 理论依据：
Romberg 算法的计算步骤如下：
（1）先求出按梯形公式所得的积分值
T 1(0)=
b -a
[f (a ) +f (b )] 2
3
（2）把区间2等分，求出两个小梯形面积之和，记为T 1(1)，即
T 1(1)=
b -a b +a
[f (a ) +f (b ) +2f ()] 42
这样由外推法可得T 2(0), T 2(0)
1
T 1(1)-() 2T 1(0)
4T 1(1)-T 1(0)。 ==
24-11-() 2
（3）把区间再等分（即22等分），得复化梯形公式T 1(2)，由T 1(1) 与T 1(2)外推可得T
(1)2
4T 1(2)-T 1(1)42T 2(1)-T 2(0)(0)=，T 3=，如此，若已算出2k 等分的复化梯形公式2
4-14-1
T 1(k ) ，则由Richardson 外推法，构造新序列 T
(k -1)
m +1
(k ) (k -1) 4m T m -T m
， m=1,2,…,l, k=1,2,…,l-m+1, =
4m -1
最后求得T l (0)+1。
(l +1) (0)(0)
（4）T l (0)≈T l (0)或
般可用如下算法：
⎧(0)b -a
[f (a ) +f (b )]⎪T 1=2⎪
2l -1⎪1b -a b -a ⎪(l ) (l -1)
T ={T +f [a +(2i -1) ]}⎨11l -1∑l
22i =12⎪
(k ) (k -1) ⎪4m T m -T m (k -1)
⎪T m +1=, m =1, 2, , l , k =1, 2, , l -m +1m
4-1⎪⎩
其具体流程如下，并全部存入第一列
(1)T
1(0)
(3)T 2(6)T 3(0) (2)T 1(1)(5)T 2(4)T 1(2)通常计算时，用固定l=N来计算，一般l=4或5即能达到要求。
3.2 C语言程序代码：
#include #include
double f(double x) //计算f(x)的值 {double z;
z=pow(3,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return(z);} main()
{ double t[20][20],s,e=0.00001,a=1,b=3; int i,j,l,k;
t[0][1]=(b-a)*(f(b)+f(a))/2; //下为romberg 算法
t[1][1]=(b-a)*(f(b)+2*f((b+a)/2)+f(a))/4; t[0][2]=(a*t[1][1]-t[0][1])/(4-1);j=3; for(l=2;fabs(t[0][j-1]-t[0][j-2])>=e;l++) {for(k=1,s=0;k
s+=f(a+(2*k-1)*(b-a)/pow(2,l));//判断前后两次所得的T(0)的差是否符合要求，如果符合精度要求则停止循环
t[l][1]=(t[l-1][1]+(b-a)*s/pow(2,l-1))/2; for(i=l-1,j=2;i>=0;i--,j++)
t[i][j]=(pow(4,j-1)*t[i+1][j-1]-t[i][j-1])/(pow(4,j-1)-1);} if(t[0][1]
printf("t=%0.6f\n",t[0][1]); else
printf("用Romberg 算法计算函数所得近似结果为：\nf(x)=%0.6f\n",t[0][j-1]);}
3.3 运行结果：
3.4 MATLAB 上机程序
function [T,n]=mromb(f,a,b,eps) if nargin
R(1,1)=(h/2)*(feval(f,a)+feval(f,b)); n=1;J=0;err=1; while (err>eps) J=J+1;h=h/2;S=0; for i=1:n x=a+h*(2*i-1); S=S+feval(f,x); end
R(J+1,1)=R(J,1)/2+h*S; for k=1:J
R(J+1,k+1)=(4^k*R(J+1,k)-R(J,k))/(4^k-1); end
err=abs(R(J+1,J+1)-R(J+1,J)); n=2*n; end R;
T=R(J+1,J+1);
format long
f=@(x)(3.^x)*(x.^1.4)*(5*x+7)*sin(x*x); [T,n]=mromb(f,1,3,1.e-5)
3.5 问题讨论：
1.Romberge 算法的优点是:把积分化为代数运算, 而实际上只需求T 1(i), 以后用递推可得. 算法简单且收敛速度快, 一般4或5次即能达到要求。
2.Romberge 算法的缺点是:对函数的光滑性要求较高，计算新分点的值时，这些数值的个数成倍增加。
3. 该程序较为复杂，涉及函数定义，有循环，而且循环中又有判断，编写时需要注意该判断条件是处于循环中，当达到要求时跳出循环，终止运算。
4. 函数的定义可放在主函数前也可在主程序后面。本程序采用的后置方式。
4. 用定步长四阶Runge-Kutta 求解
⎧dy 1/dt =1⎪dy /dt =y
3⎪2
⎪⎪dy 3/dt =1000-1000y 2-100y 3
⎨
y (0) =0⎪1
⎪y 2(0) =0⎪⎪⎩y 3(0) =0
h =0.0005，打印y i (0.025) ， y i (0.045) ， y i (0.085) ， y i (0.1) ，(i =1，2，3) 4.1 理论依据：
Runge_Kutta采用高阶单步法，这里不是先按Taylor 公式展开，而是先写成t n 处附近的值的线性组合（有待定常数）再按Taylor 公式展开，然后确定待定常数，这就是Runge-Kutta 法的思想方法。
本题采用四阶古典的Runge-Kutta 公式：
Y n +1=Y n +[K 1+3K 2+3K 3+K 4]/8K 1=hF (x n , Y n )
K 2=hF (x n +h /3, Y n +hK 1/3) K 3=hF (x n +2h /3, Y n +hK 1/3+hK 2) K 4=hF (x n +h , Y n +hK 1-hK 2+hK 3)
4.2 C语言程序代码：
#include
void fun(double x[4],double y[4],double h) {y[1]=1*h; y[2]=x[3]*h;
y[3]=(1000-1000*x[2]-100*x[2]-100*x[3])*h; //微分方程向量函数} void main()
{ double Y[5][4],K[5][4],m,z[4],e=0.0005; double y[5]={0,0.025,0.045,0.085,0.1}; int i,j,k;
for(i=1;i
for(j=1;j
for(k=1;k
{for(m=y[k-1];m
z[i]=Y[k][i]+e*K[2][i]/2; //依此求K1,K2K3的值 fun(z,K[2],e); for(i=1;i
z[i]=Y[k][i]+e*K[2][i]/2; fun(z,K[3],e); for(i=1;i
z[i]=Y[k][i]+e*K[3][i]; fun(z,K[4],e); for(i=1;i
Y[k][i]=Y[k][i]+(K[1][i]+2*K[2][i]+2*K[3][i]+K[4][i])/6; // 求Yi[N+1]的值} if(k!=5)
for(i=1;i
Y[k+1][i]=Y[k][i];} printf("计算结果:\n"); for(i=1;i
{printf("y%d[%4.3f]=%-10.8f,",j,y[i],Y[i][j]); if(j==3) printf("\n");} printf("\n");} }
4.3 运行结果：
4.4 问题讨论：
1. 定步长四阶Runge-kutta 方法是一种高阶单步法法稳定，精度较高，误差小且程序相对简单，存储量少。不必求出起始点的函数值，可根据精度的要求修改步长，不会由于起始点的误差造成病态。
2. 本程序可以通过修改主程序所调用的函数中的表达式来实现对其它函数的任意初值条件求微分计算。
3. 程序中运用了大量的for 循环语句，因为该公式中涉及大量的求和，且有不同的函数和对不同的数值求值，编程稍显繁琐。所以编写过程中一定要注意各循环的次数，以免出错。
5.
2.115237 -1.061074 1.112336 -0.113584 0.718719 1.742382 3.067813 -2.031743⎤⎡12.38412
⎢2.115237⎥ 19.141823 -3.125432 -1.012345 2.189736 1.563849 -0.784165 1.112348 3.123124⎢⎥⎢-1.061074 -3.125432 15.567914 3.123848 2.031454 1.836742 -1.056781 0.336993 -1.010103⎥⎢⎥1.112336 -1.012345 3.123848 27.108437 4.101011 -3.741856 2.101023 -0.71828 -0.037585⎢⎥
⎥A =⎢-0.113584 2.189736 2.031454 4.101011 19.897918 0.431637 -3.111223 2.121314 1.784317
⎢⎥0.718719 1.563849 1.836742 -3.741856 0.431637 9.789365 -0.103458 -1.103456 0.238417⎢⎥⎢1.742382⎥ -0.784165 -1.056781 2.101023 -3.111223 -0.103458 14.7138465 3.123789 -2.213474⎢⎥
1.112348 0.336993 -0.71828 2.121314 -1.103456 3.123789 30.719334 4.446782⎢3.067813⎥
⎢-2.031743⎥ 3.123124 -1.010103 -0.037585 1.784317 0.238417 -2.213474 4.446782 40.00001⎣⎦
b =(2.1874369 33.992318 -25.173417 0.84671695 1.784317 -86.612343 1.1101230 4.719345 -5.6784392) T
用列主元消去法求解Ax=b。
5.1 理论依据：
列主元素消元法是在应用Gauss 消元法的基础上，凭借长期经验积累提出的，是线性方程组一般解法，目的是为避免在消元计算中使误差的扩大，甚至严重损失了有效数字使数据失真，而在每次初等变换前对矩阵作恰当的调整，以提高Gauss 消元法的数字稳定性，进而提高计算所得数据的精确度。即在每主列中取绝对值最大的元素作主元，再做对应的行交换然后消元求解的办法。具体做法如下：
将方阵A 和向量b 写成C=（A ，b ）。将C 的第1列中第1行的元素与其下面的此列的元素逐一进行比较，找到最大的元素c j 1，将第j 行的元素与第1行的元素进行交换，然后通过行变换，将第1列中第2到第n 个元素都消成0。将变换后的矩阵C (1)的第二列中第二行的元素与其下面的此列的元素逐一进行比较，找到最大的元素c k (1)2，将第k 行的元素与第2行的元素进行交换，然后通过行变换，将第2列中第3到第n 个元素都消成0。以此方法将矩阵的左下部分全都消
成0后再求解。最终形式如下：
(n ) ⎛a 11 0
(A,b）～
0⎝
*
(n )
g 1⎫a 1n
⎪ ⎪ ⎪
⎪(n )
a nn g n ⎭
5.2 C语言程序代码
（1）比较该列的元素的绝对值的大小，将绝对值最大的元素通过行变换使其位于主对角线上；
（2）进行高斯消去法变换，把系数矩阵化成上三角形，然后回代求
#include "math.h" #include "stdio.h"
void Householder(double A[9][9]);
void expunction(double A[9][9],double b[9],double x[9]); void main() {double A[9][9]={
{12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743},
{2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3.123124},
{-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1.010103},
{1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585},
{-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317},
{0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417},
{1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713847,3.123789,-2.213474},
{3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782},
{-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40.00001}}; double b[9]=
{2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392}; double x[9]={0.0}; int i,j;
Householder(A);
printf("\n The Results of X are:\n"); expunction(A,b,x); for(i=1;i
printf("X%1d=%f\n",i,x[i-1]);}
void Householder(double A[9][9]) {double q[9],u[9],y[9],s,a,kr; int i,j,k; for(i=0;i
for(j=i+1;j
a=s*s+fabs(A[i+1][i])*s; for(j=0;j
else if(j==i+1) u[j]=A[j][i]+A[j][i]/fabs(A[j][i])*s; else if(j>i+1) u[j]=A[j][i];} for(k=0;k
{y[k]=0; for(j=0;j
for(k=0;k
for(k=0;k
A[k][j]-=u[k]*q[j]+u[j]*q[k];} }
}
void expunction(double A[9][9],double b[9],double x[9]) {int i,j,k; double B[9][10]; double z[3];
double t1=0,t2=0,t3=0; for(i=0;iA[i][i]) {for(j=i,k=0;j
z[k]=A[i][j];A[i][j]=A[i+1][j];A[i+1][j]=z[k]; t1=b[i];b[i]=b[i+1];b[i+1]=t1;} t2=A[i+1][i]; for(j=i;j
A[i+1][j]=A[i+1][j]-A[i][j]*t2/A[i][i]; b[i+1]=b[i+1]-b[i]*t2/A[i][i];} x[8]=b[8]/A[8][8]; for(i=7;i>=0;i--)
{for(j=i+1;j
5.3 运行结果
5.4 MATLAB 上机程序
unction [x]=mgauss2(A,b,flag) if nargin
[ap,p]=max(abs(A(k:n,k))); p=p+k-1; if p>k
A([k p],:)=A([p k],:); b([k p],:)=b([p k],:); end
m=A(k+1:n,k)/A(k,k);
A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n); b(k+1:n)=b(k+1:n)-m*b(k); A(k+1:n,k)=zeros(n-k,1); if flag~=0,Ab=[A,b],end end
x=zeros(n,1); x(n)=b(n)/A(n,n);
for k=n-1:-1:1
x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);
end
format long
A=[12.38412,2.115237,-1.061074,1.112336,-0.113584,0.718719,1.742382,3.067813,-2.031743;
2.115237,19.141823,-3.125432,-1.012345,2.189736,1.563849,-0.784165,1.112348,3. 123124;
-1.061074,-3.125432,15.567914,3.123848,2.031454,1.836742,-1.056781,0.336993,-1. 010103;
1.112336,-1.012345,3.123848,27.108437,4.101011,-3.741856,2.101023,-0.71828,-0.037585;
-0.113584,2.189736,2.031454,4.101011,19.897918,0.431637,-3.111223,2.121314,1.784317;
0.718719,1.563849,1.836742,-3.741856,0.431637,9.789365,-0.103458,-1.103456,0.238417;
1.742382,-0.784165,-1.056781,2.101023,-3.111223,-0.103458,14.713846,3.123789,-2.213474;
3.067813,1.112348,0.336993,-0.71828,2.121314,-1.103456,3.123789,30.719334,4.446782;
-2.031743,3.123124,-1.010103,-0.037585,1.784317,0.238417,-2.213474,4.446782,40. 00001];
b=[2.1874369,33.992318,-25.173417,0.84671695,1.784317,-86.612343,1.1101230,4.719345,-5.6784392]';
x=mgauss2(A,b);x=x'
5.5问题讨论：
1. 输入矩阵可用循环比较指令确认输入的是否为对称矩阵。
2. 循环体中的累加值注意初始化零。
3. 注意公式中的下标从一开始，数组中的下标从零开始。
4. 在编程中数组的角标从0开始与数学中的起始脚表不同，编程时必须注意。
5. 在给数组元素通过表达式赋值时要注意原始赋值的覆盖问题。
20