实验3非线性方程AX0解法.docx
《实验3非线性方程AX0解法.docx》由会员分享，可在线阅读，更多相关《实验3非线性方程AX0解法.docx（33页珍藏版）》请在冰豆网上搜索。
实验3非线性方程AX0解法
线性方程组AX=B的数值解法
1.实验描述
1．，2,3:
通过矩阵可表示立方体的坐标位置，与另一矩阵相乘可实现立方体坐标位置进行变换
2．:
不通过行变换就可以解决三角线性方程。
：
将单位矩阵表示成列矩阵，通过对目标矩阵别离求解得出列矩阵从而取得目标矩阵的逆矩阵。
4.：
将单位矩阵表示成列矩阵，通过度解成上下三角矩阵对目标矩阵别离求解得出列矩阵从而取得目标矩阵的逆矩阵。
应用程序求解基尔霍夫电流。
应用高斯-赛德尔迭代法求解带状方程。
2.实验内容
.
单位立方体位于第一卦限，一个极点在原点。
第一，以角度
沿y轴旋转立方体，然后再以角度
沿z轴旋转立方体。
求旋转后立方体的8个极点的坐标，并与例的结果比较。
它们的区别是什么？
试通过矩阵一般不知足互换律的事实进行解释。
利用plot3命令画出3个图形。
.
设单位立方体位于第一卦限，其中一个极点位于坐标原点。
第一以角度
沿x轴旋转立方体，然后再以角度
沿z轴旋转立方体。
求旋转后立方体的8个极点的坐标。
利用plot3画出这3个立方体。
.
四面体的坐标为（0,0,0），（1,0,0），（0,1,0），（0,0,1）。
第一以弧度沿y轴旋转，然后再以弧度沿z轴旋转，最后以弧度沿x轴旋转。
求旋转后的极点坐标。
利用plot3画出这4个立方体。
.许多科学应用包括的矩阵带有很多零。
在实际情形中很重要的三角形线性方程组有如下形式：
……
构造一个程序求解三角形线性方程组。
可假定不需要变换。
而且可用第k行消去第k+1行的
。
下面的习题虽然是针对3x3维矩阵的，但其概念可用于NxN维矩阵。
若是矩阵A非奇异，则
存在。
而且
。
设
，
，
是
的列，而
，
，
是I的列。
方程
可表示为
则上式等价于三个线性方程组
，
，
如此求
等价于求解三个线性方程组。
利用程序或上题中的程序求解下面每一个矩阵的逆。
通过计算
和利用命令inv（A）检查答案，并解释可能的不同。
（a）
（b）
修改程序，使得它能够通过重复求解N个线性方程组
其中J=1，2，…，N
来取得
，
则
、
而且
，保证对LU分解只计算一次。
基尔霍夫电压定律说明电路网络中任意单向闭路的电压和为零。
现有如下电线路性方程：
若是
利用程序求解电流
利用高斯-赛德尔迭代法求解下列带状方程。
3.实验结果及分析
:
算法:
（1）令X=zeros（8,3）;X（[5:
8,11,12,15,16,18,20,22,24]）=1;d=[1243156875624873];i=0。
（2）判断i>100是不是成立，若成立，执行步骤（3）；若不成立，r1=[cos（i*pi/600）-sin（i*pi/600）0;010;-sin（i*pi/600）0cos（i*pi/600）];U=X*r1';plot3（U（d,1）,U（d,2）,U（d,3））;drawnow，i=i+1,返回步骤
（2）.\
（3）i=0.
（4）判断i>100是不是成立，若成立，执行步骤（4）；若不成立，r2=[cos（i*pi/400）-sin（i*pi/400）0;sin（i*pi/400）cos（i*pi/400）0;001];W=U*r2';plot3（W（d,1）,W（d,2）,W（d,3））;drawnow，i=i+1,返回步骤（3）.
（5）subplot（2,2,1）;plot3（X（d,1）,X（d,2）,X（d,3））subplot（2,2,2）;plot3（U（d,1）,U（d,2）,U（d,3））;subplot（2,2,3）；
plot3（W（d,1）,W（d,2）,W（d,3））;xlabel（'x'）；ylabel（'y'）；zlabel（'z'）；view（3）;rotate3d。
图（左上）初始立方体
图（右上）
，沿y轴旋转
图（左下）
沿z轴旋转
表：
第一次旋转后立方体的坐标
X
0
0
Y
0
0
0
0
Z
0
0
表：
第二次旋转后立方体的坐标
X
0
Y
0
Z
0
0
算法:
（1）令X=zeros（8,3）;X（[5:
8,11,12,15,16,18,20,22,24]）=1;d=[1243156875624873];i=0。
（2）判断i>100是不是成立，若成立，执行步骤（3）；若不成立，r1=[100;0cos（i*pi/1200）-sin（i*pi/1200）;0sin（i*pi/1200）cos（i*pi/1200）];;U=X*r1';plot3（U（d,1）,U（d,2）,U（d,3））;drawnow，i=i+1,返回步骤
（2）.\
（3）i=0.
（4）判断i>100是不是成立，若成立，执行步骤（4）；若不成立，r2=[cos（i*pi/600）-sin（i*pi/600）0;sin（i*pi/600）cos（i*pi/600）0;001];plot3（W（d,1）,W（d,2）,W（d,3））;drawnow，i=i+1,返回步骤（3）.
（5）subplot（2,2,1）;plot3（X（d,1）,X（d,2）,X（d,3））subplot（2,2,2）;plot3（U（d,1）,U（d,2）,U（d,3））;subplot（2,2,3）；
plot3（W（d,1）,W（d,2）,W（d,3））;xlabel（'x'）；ylabel（'y'）；zlabel（'z'）；view（3）;rotate3d。
图（左上）单位立方体
图（右上）第一次旋转
图（左下）第2次旋转
表：
第一次旋转后立方体坐标
X
0
1
0
0
1
1
0
1
Y
0
0
Z
0
0
表：
第二次旋转后立方体坐标
X
0
Y
0
Z
0
0
算法：
（1）令X=zeros（8,3）;X（[5:
8,11,12,15,16,18,20,22,24]）=1;d=[1243156875624873];i=0。
（2）判断i>100是不是成立，若成立，执行步骤（3）；若不成立，r1=[100;0cos（i*pi/1200）-sin（i*pi/1200）;0sin（i*pi/1200）cos（i*pi/1200）];;U=X*r1';plot3（U（d,1）,U（d,2）,U（d,3））;drawnow，i=i+1,返回步骤
（2）.\
（3）i=0.
（4）判断i>100是不是成立，若成立，执行步骤（4）；若不成立，r2=[cos（i*pi/600）-sin（i*pi/600）0;sin（i*pi/600）cos（i*pi/600）0;001];W=U*r2'plot3（W（d,1）,W（d,2）,W（d,3））;drawnow，i=i+1,返回步骤（3）.
（5）i=0.
（6）判断i>100是不是成立，若成立，执行步骤（7）；若不成立，r3=[100;0cos（i*100）-sin（i*100）;0sin（i*100）cos（i*100）];T=W*r3';
plot3（T（d,1）,T（d,2）,T（d,3））;drawnow，i=i+1,返回步骤（6）
（7）subplot（2,2,1）;plot3（X（d,1）,X（d,2）,X（d,3））subplot（2,2,2）;plot3（U（d,1）,U（d,2）,U（d,3））;subplot（2,2,3）；
plot3（W（d,1）,W（d,2）,W（d,3））;subplot（2,2,4）;plot3（T（d,1）,T（d,2）,T（d,3））;xlabel（'x'）；ylabel（'y'）；zlabel（'z'）；view（3）;rotate3d。
图（左上）正四面体
图（右上）第一次旋转
图（左下）第二次旋转
图（右下）第4次旋转
表四面体坐标
X
0
1
0
0
Y
0
0
1
0
Z
0
0
0
1
表第一次旋转后坐标
X
0
0
Y
0
0
0
Z
0
0
表第二次旋转后坐标
X
0
Y
0
Z
0
0
表第三次旋转后坐标
X
0
Y
0
Z
0
算法：
（1）输入A,B,P,delta,max1，令N=length（B）;k=1.
（2）判断k>max1是不是成立，若成立，输出结果；若不成立，执行步骤（3）。
（3）令j=1,判断j>N是不是成立，若成立执行步骤（6）；若不成立，执行步骤（4）。
（4）判断j==1是不是成立，若成立，X
（1）=（B
（1）-A（1,2）*P
（2））/A（1,1），j=j+1,
执行步骤（3）;若不成立，v执行步骤（5）。
（5）判断j==N是不是成立，若成立，X（N）=（B（N）-A（N,N-1）*（X（N-1））'）/A（N,N），j=j+1,
执行步骤（3）;若不成立，X（j）=（B（j）-A（j,j-1）*X（j-1）'-A（j,j+1）*P（j+1））/A（j,j），j=j+1,执行步骤（3）。
（6）令err=abs（norm（X'-P））;relerr=err/（norm（X）+eps）;P=X';
（7）判断（err（2）.
流程图：
算法：
（1）输入A,B，令[NN]=size（A）;X=zeros（N,1）;C=zeros（1,N+1）;Aug=[AB];p=1。
（2）判断p>N-1是不是成立，若不成立，输出结果；若成立，[Y,j]=max（abs（Aug（p:
N,p）））;
C=Aug（p,:
）;Aug（p,:
）=Aug（j+p-1,:
）;Aug（j+p-1,:
）=C;
（3）判断Aug（p,p）==0是不是成立，若成立，输出结果；若不成立，执行步骤（4）。
（4）令k=p+1，判断k>N是不是成立，若成立，p=p+1,返回步骤
（2）;若不成立，m=Aug（k,p）/Aug（p,p）;
Aug（k,p:
N+1）=Aug（k,p:
N+1）-m*Aug（p,p:
N+1），k=k+1,返回步骤（4）。
（5）X=backsub（Aug（1:
N,1:
N）,Aug（1:
N,N+1）），输出结果。
a）
（b）
修改程序主要变更为列矩阵B变成了单位矩阵E
算法：
（1）输入A,E；令[N,N]=size（A）;X=zeros（N,1）;Y=zeros（N,1）;C=zeros（1,N）;t=1:
N;p=1。
（2）判断p=1>N-1是不是成立，若成立，执行步骤（6）；若不成立，执行步骤（3）。
（3）[max1,j]=max（abs（A（p:
N,p）））;C=A（p,:
）;A（p,:
）=A（j+p-1,:
）;A（j+p-1,:
）=C;d=t（p）;t（p）=t（j+p-1）;t（j+p-1）=d;
（4）判断A（p,p）==0是不是成立，若成立，执行步骤（12）；若不成立，k=p+1，执行步骤（5）。
（5）判断k>N是不是成立，若成立，返回到步骤
（2）；若不成立，mult=A（k,p）/A（p,p）;
A（k,p）=mult;A（k,p+1:
N）=A（k,p+1:
N）-mult*A（p,p+1:
N）;k=k+1;返回到步骤（5）。
（6）J=1。
（7）判断J>N是不是成立，若成立，执行步骤（12）；若不成立，k=2,Y
（1）=E（t
（1）,J）;执行步骤（8）
（8）判断k>N是不是成立，若成立，执行步骤（9）；若不成立，Y（k）=E（t（k）,J）-A（k,1:
k-1）*Y（1:
k-1），k=k+1,返回步骤（8）。
（9）X（N）=Y（N）/A（N,N）;k=N-1.
（10）判断k>1是不是成立，若成立，执行步骤（11）；若不成立，X（k）=（Y（k）-A（k,k+1:
N）*X（k+1:
N））/A（k,k）;
k=k-1;执行步骤（10）.
（11）m（1:
N,J）=X（1:
N）;返回步骤（7）.
（12）输出结果。
流程图：
算法：
（1）输入A,E；令[N,N]=size（A）;X=zeros（N,1）;Y=zeros（N,1）;C=zeros（1,N）;t=1:
N;p=1。
（2）判断p=1>N-1是不是成立，若成立，执行步骤（6）；若不成立，执行步骤（3）。
（3）[max1,j]=max（abs（A（p:
N,p）））;C=A（p,:
）;A（p,:
）=A（j+p-1,:
）;A（j+p-1,:
）=C;d=t（p）;t（p）=t（j+p-1）;t（j+p-1）=d;
（4）判断A（p,p）==0是不是成立，若成立，执行步骤（11）；若不成立，k=p+1，执行步骤（5）。
（5）判断k>N是不是成立，若成立，返回到步骤
（2）；若不成立，mult=A（k,p）/A（p,p）;
A（k,p）=mult;A（k,p+1:
N）=A（k,p+1:
N）-mult*A（p,p+1:
N）;k=k+1;返回到步骤（5）。
（6）k=2,Y
（1）=E（t
（1）,J）;执行步骤（7）
（7）判断k>N是不是成立，若成立，执行步骤（8）；若不成立，Y（k）=E（t（k）,J）-A（k,1:
k-1）*Y（1:
k-1），k=k+1,返回步骤（7）。
（8）X（N）=Y（N）/A（N,N）;k=N-1.
（9）判断k>1是不是成立，若成立，执行步骤（10）；若不成立，X（k）=（Y（k）-A（k,k+1:
N）*X（k+1:
N））/A（k,k）;
k=k-1;执行步骤（9）.
（10）输出结果。
流程图：
表3电流
（a）
3
5
1
（b）
8
1
（c）
4
3
-1
.
算法：
（1）输入A,B,P,delta,max1，令N=length（B）;k=1.
（2）判断k>max1是不是成立，若成立，输出结果；若不成立，执行步骤（3）。
（3）令j=1,判断j>N是不是成立，若成立执行步骤（8）；若不成立，执行步骤（4）。
（4）判断j==1是不是成立，若成立，X
（1）=（B
（1）-A（1,2）*P
（2））/A（1,1），j=j+1,
执行步骤（3）;若不成立，执行步骤（5）。
（5）判断j==2是不是成立，若成立，X
（2）=（B
（2）-A（2,1）*X
（1）'-A（2,3:
4）*P（3:
4））/A（2,2）;
执行步骤（3）;若不成立，执行步骤（6）。
（6）判断j==N-1是不是成立，若成立，X（N-1）=（B（N-1）-A（N-1,N-3:
N-2）*X（N-3:
N-2）'-A（N-1,N）*P（N））/A（N-1,N-1）;执行步骤（3）;若不成立，执行步骤（7）。
（7）判断j==N是不是成立，若成立，X（N）=（B（N）-A（N,N-1）*（X（N-1））'）/A（N,N），j=j+1,
执行步骤（3）;若不成立，X（j）=（B（j）-A（j,j-1）*X（j-1）'-A（j,j+1）*P（j+1））/A（j,j），j=j+1,执行步骤（3）。
（8）令err=abs（norm（X'-P））;relerr=err/（norm（X）+eps）;P=X';
（9）判断（err执行步骤
（2）.
表4.通太高斯赛德尔迭代取得的x的解
1~10
0.
0.
0.
0.
0.
0.
11~20
0.
0.
0.
0.
21~30
0.
0.
31~40
0.
0.
0.
41~50
0.
0.
0.
0.
0.
0.
0.
4.结论
2,3
1）利用矩阵乘的方式可对已知三维图形进行坐标变换，实现该图形在空间的旋转。
2）矩阵乘一般不知足互换律，由于利用矩阵乘的方式对三维图形进行坐标变换，例如V=R1*U;W=R2*V，现在
;故在进行多次旋转时，需要分步进行，不可一步完成旋转。
7
通过对系数矩阵的增广矩阵进行高斯消元和回带容易患到线性方程组的解，同时，利用这种方式能够求得矩阵的逆。
4
若是系数矩阵能分解为LU的形式，其中L为下三角矩阵，U为上三角矩阵，通过对系数矩阵的分解再求解可应用简单的迭代进行求解x。
.
求解带状线性方程组的解可利用高斯-赛德尔迭代法。
附件（代码）：
:
X=zeros（8,3）;
X（[5:
8,11,12,15,16,18,20,22,24]）=1;
d=[1243156875624873];
fori=0:
1:
100
r1=[cos（i*pi/600）-sin（i*pi/600）0;010;-sin（i*pi/600）0cos（i*pi/600）];
U=X*r1';
plot3（U（d,1）,U（d,2）,U（d,3））;
drawnow
end
fori=0:
1:
100
r2=[cos（i*pi/400）-sin（i*pi/400）0;sin（i*pi/400）cos（i*pi/400）0;001];
W=U*r2';
plot3（W（d,1）,W（d,2）,W（d,3））;
drawnow
end
subplot（2,2,1）;
plot3（X（d,1）,X（d,2）,X（d,3））;
subplot（2,2,2）;
plot3（U（d,1）,U（d,2）,U（d,3））;
subplot（2,2,3）;
plot3（W（d,1）,W（d,2）,W（d,3））;
xlabel（'x'）
ylabel（'y'）
zlabel（'z'）
view（3）;rotate3d;
X=zeros（8,3）;
X（[5:
8,11,12,15,16,18,20,22,24]）=1;
d=[1243156875624873];
fori=0:
1:
100
r1=[100;0cos（i*pi/1200）-sin（i*pi/1200）;0sin（i*pi/1200）cos（i*pi/1200）];
U=X*r1';
plot3（U（d,1）,U（d,2）,U（d,3））;
drawnow
end
fori=0:
1:
100
r2=[cos（i*pi/600）-sin（i*pi/600）0;sin（i*pi/600）cos（i*pi/600）0;001];
W=U*r2';
plot3（W（d,1）,W（d,2）,W（d,3））;
drawnow
xlabel（'x'）
ylabel（'y'）
zlabel（'z'）
end
subplot（2,2,1）;
plot3（X（d,1）,X（d,2）,X（d,3））;
subplot（2,2,2）;
plot3（U（d,1）,U（d,2）,U（d,3））;
subplot（2,2,3）;
plot3（W（d,1）,W（d,2）,W（d,3））;
view（3）;rotate3d;
xlabel（'x'）
ylabel（'y'）
zlabel（'z'）
X=zeros（4,3）;
X（[2,7,12]）=1;
d=[12341324];
fori=0:
1:
100
r1=[cos（i*100）-sin（i*100）0;010;-sin（i*100）0cos（i*100）];
U=X*r1';
plot3（U（d,1）,U（d,2）,U（d,3））;
drawnow
end
fori=0:
1:
100
r2=[cos（-i*100）-sin（-i*100）0;sin（-i*100）cos（-i*100）0;001];
W=U*r2';
plot3（W（d,1）,W（d,2）,W（d,3））;
drawnow
end
fori=0:
1:
100
r3=[100;0cos（i*100）-sin（i*100）;0sin（i*100）cos（i*100）];
T=W*r3';
plot3（T（d,1）,T（d,2）,T（d,3））;
drawnow
end
subplot（2,2,1）;
plot3（X（d,1）,X（d,2）,X（d,3））;
subplot（2,2,2）;
plot3（U（d,1）,U（d,2）,U（d,3））;
subplot（2,2,3）;
plot3（W（d,1）,W（d,2）,W（d,3））;
subplot（2,2,4）;
plot3（T（d,1）,T（d,2）,T（d,3））;
view（3）;rotate3d;
xlabel（'x'）
ylabel（'y'）
zlabel（'z'）
:
functionX=gseid2（A,B,P,delta,max1）
%input--AisanN*Nnonsingularmatrix
%-BisanN*1matrix
%-PisanN*1matrix
%-deltaisthetoleranceforP
%-max1isthemaximumnumberofiterations
%output-XisanN*1matrix:
thegauss-seidelapproximationtothesolution
%ofAX=B
N=length（B）;
fork=1:
max1
forj=1:
N
ifj==1
X
（1）=
升级会员 