实验3 - 非线性方程AX=0的解法

更新时间：2023-10-30 12:20:01 阅读量：综合文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

实验3r原则推荐度：
相关推荐

《数值计算方法》实验报告

线性方程组AX=B的数值解法

1.实验描述

1．P93.1，2,3:通过矩阵可表示立方体的坐标位置，与另一矩阵相乘可实现立方体坐标位置进行变换

2．P108.1:不通过行变换就能解决三角线性方程。

3.p109.7：将单位矩阵表示成列矩阵，通过对目标矩阵分别求解得出列矩阵从而得到目标矩

阵的逆矩阵。

4. 120.3：将单位矩阵表示成列矩阵，通过分解成上下三角矩阵对目标矩阵分别求解得出列矩阵从而得到目标矩阵的逆矩阵。 5.p120.4:应用程序3.3求解基尔霍夫电流。 6.p129.4:应用高斯-赛德尔迭代法求解带状方程。

2.实验内容

P93.1.

单位立方体位于第一卦限，一个顶点在原点。首先，以角度再以角度

?沿y轴旋转立方体，然后6?沿z轴旋转立方体。求旋转后立方体的8个顶点的坐标，并与例3.10的结果比较。4它们的区别是什么？试通过矩阵一般不满足交换律的事实进行解释。使用plot3命令画出3个图形。 P93.2.

?设单位立方体位于第一卦限，其中一个顶点位于坐标原点。首先以角度沿x轴旋转立

12?方体，然后再以角度沿z轴旋转立方体。求旋转后立方体的8个顶点的坐标。使用plot3

6画出这3个立方体。 P93.3.

四面体的坐标为（0,0,0），（1,0,0），（0,1,0），（0,0,1）。首先以弧度0.15沿y轴旋转，然后再以弧度-1.5沿z轴旋转，最后以弧度2.7沿x轴旋转。求旋转后的顶点坐标。使用plot3画出这4个立方体。 P108.1

.许多科学应用包含的矩阵带有很多零。在实际情况中很重要的三角形线性方程组有如下形式：

d1x1?c1x? b12《数值计算方法》实验报告 2

a1x1?d2x2?c2x?3 b a2x2?d3x3?c3x?4 b ……

aN?2xN?2?dN?1xN?1?cN?1xN? aN?1xN?1?dNxN?bN

构造一个程序求解三角形线性方程组。可假定不需要变换。而且可用第k行消去第k+1行的xk。

p109.7

下面的习题虽然是针对3x3维矩阵的，但其概念可用于NxN维矩阵。如果矩阵A非奇异，则A?1存在。而且AA?1?I。设C1，而E1，方程AA?1?IC2，C3是A?1的列，E2，E3是I的列。可表示为A?C1,C2,C3???E1,E2,E3?

则上式等价于三个线性方程组AC1?E1，AC2?E2，AC3?E3 这样求A?1等价于求解三个线性方程组。

使用程序2.2或上题中的程序求解下面每个矩阵的逆。通过计算 AA?1 和使用命令inv(A)检查答案，并解释可能的差异。

N?b?120240?140??16?201???1201200?27001680???? （a）?325? （b）??240?27006480?4200??????1?10?42002800???1401680120.3

修改程序3.3，使得它可以通过重复求解N个线性方程组

ACJ?EJ 其中J=1，2，…，N 来得到A?1，

则 A?C1C2...CN???E1E2...EN?、

而且 A?1??C1C2...CN?，保证对LU分解只计算一次。 p120.4

基尔霍夫电压定律说明电路网络中任意单向闭路的电压和为零。现有如下电路线性方程：

《数值计算方法》实验报告

(R1?R3?R4)I1?R3I2?R4I3?E1R3I1?(R2?R3?R5)I2?R5I3?E2 R4I1?R5I2?(R4?R5?R6)I3?0 如果

(a)R1?1,R2?1,R3?2,R4?1,R5?2,R6?4,E1?23,E2?29(b)R1?1,R2?0.75,R3?1,R4?2，R5?1,R6?4,E1?12,E2?21.5(c)R1?1,R2?2,R3?4,R4?3,R5?1,R6?5,E1?41,E2?38 使用程序求解电流I1,I2,I3。 p129.4

利用高斯-赛德尔迭代法求解下列带状方程。

12x1?2x2?x3?5?2x1?12x2?2x3?x4?5x1?2x2?12x3?2x4?x5?5x2?2x3?12x4?2x5?x6?5??????x46?2x47?12x48?2x49?x50?5x47?2x48?12x49?2x50?5x48?2x49?12x50?5

3.实验结果及分析

P93.1: 算法:

(1) 令X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0。 (2) 判断i>100是否成立，若成立，执行步骤（3）；若不成立， r1=[cos(i*pi/600) -sin(i*pi/600) 0;0 1 0;-sin(i*pi/600) 0 cos(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;0 0 1];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。

《数值计算方法》实验报告 4

start X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0 i=i+1 Y i>100 N r1=[cos(i*pi/600),-sin(i*pi/600),0;0,1,0;-sin(i*pi/600,0 ,cos(i*pi/600)];U=X*r1';plot3(U(d,1),U(d,2),U(d,3));drawnow i=0. i=i+1 Y i>100 N r2=[cos(i*pi/400),-sin(i*pi/400),0;sin(i*pi/400),cos(i*pi/400),0;0,0,1];W=U*r2';plot3(W(d,1),W(d,2),W(d,3));drawnow 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 output end

《数值计算方法》实验报告

10.5010.5000.521010.50-101121020y-2-2x0z2

图1.1(左上) 初始立方体

? 图1.2（右上）V?Ry()U，沿y轴旋转

6?图1.3 （左下）W?Rx()V沿z轴旋转

4 表1.1：第一次旋转后立方体的坐标

Z 0 -0.5000 X Y 0 0 0.8660 0.8660 0 1.0000 0 0.5000 0.5000 1.0000 0 0.8660 0.8660 1.3660 1.3660 0 1.0000 1.0000 -0.5000 0 0.3660 0.3660

表1.2：第二次旋转后立方体的坐标

X 0 0.6124 -0.094-0.7071 -0.3530.3536

7 6 Y 0 0.6124 1.3195 0.7071 1.0607 0.3536

Z 0 -0.5000 -0.5000 0 0.8660 0.8660 0.9659 0.2588 0.9659 1.6730 0.3660 0.3660

《数值计算方法》实验报告 6

P93.2 算法:

(1) 令X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0。 (2) 判断i>100是否成立，若成立，执行步骤（3）；若不成立， r1=[1 0 0;0 cos(i*pi/1200)

-sin(i*pi/1200) ;0 sin(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;0 0 1];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。

《数值计算方法》实验报告

start 7

X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0 i=i+1 Y i>100 N r1=[1,0,0;0,cos(i*pi/1200),-sin(i*pi/1200);0,sin(i*pi/1200),cos(i*pi/1200)];U=X*r1';plot3(U(d,1),U(d,2),U(d,3));drawnow i=0. i=i+1 Y i>100 N r2=[cos(i*pi/600),-sin(i*pi/600),0;sin(i*pi/600),cos(i*pi/600),0;0,0,1];plot3(W(d,1),W(d,2),W(d,3));W=U*r2';plot3(W(d,1),W(d,2),W(d,3));drawnow 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 output end 《数值计算方法》实验报告 8

10.5010.5000.521010-100.51121020y-2-1x0z1

图2.1（左上）单位立方体图2.2 （右上）第一次旋转图2.3 （左下）第2次旋转

表2.1：第一次旋转后立方体坐标

表2.2：第二次旋转后立方体坐标

X Y Z 0 0 0 0.8660 0.5000 0 -0.4830 0.1294 0.3831 0.8365 0.2588 -0.2241 0.9659 0.2588 0.9659 1.2247 1.2247 1.3365 0.9954 -0.3536 0.2759 0.6124 0.5125 1.1124 X Y Z 0 0 0 1 0 0 0 0.9659 0.2588 0 1 1 0 1 0.7071 1.2247 -0.2588 0.9659 0.9659 0.2588 -0.2588 0.7071 0.9659 1.2247

《数值计算方法》实验报告

P93.3 算法：

(1) 令X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0。 (2) 判断i>100是否成立，若成立，执行步骤（3）；若不成立， r1=[1 0 0;0 cos(i*pi/1200) -sin(i*pi/1200) ;0 sin(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;0 0 1]; 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=[1 0 0;0 cos(i*2.7/100) -sin(i*2.7/100) ;0 sin(i*2.7/100) cos(i*2.7/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)); ylabel('y')；zlabel('z')；view(3); rotate3d。

xlabel('x')；

《数值计算方法》实验报告 10

start X=zeros(8,3);X([5:8,11,12,15,16,18,20,22,24])=1;d=[1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3]; i=0 i=i+1 i>100 N Y r1=[1,0,0;0,cos(i*pi/1200),-sin(i*pi/1200);0,sin(i*pi/1200),cos(i*N pi/1200)];U=X*r1';plot3(U(d,1),U(d,2),U(d,3));drawnow i=0. i=i+1 i>100 N r2=[cos(i*pi/600),-sin(i*pi/600,0;sin(i*pi/600),cos(i*pi/600),0;0 ,0,1];W=U*r2';plot3(W(d,1),W(d,2),W(d,3));drawnow N Y i=0. i=i+1 i>100 N Y r2=[cos(i*pi/600),-sin(i*pi/600,0;sin(i*pi/600),cos(i*pi/600),0;0 ,0,1];W=U*r2';plot3(W(d,1),W(d,2),W(d,3));drawnow 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 output end

《数值计算方法》实验报告

10.5010.5000.510-110.50-101110-120-200.110-110y-10x0.1z0.20.2

图1.1（左上）正四面体图1.2 （右上）第一次旋转图3.3 （左下）第二次旋转图3.4 （右下）第4次旋转 X Y Z 0 0 0 表3.1 四面体坐标 1 0 0 0 1 0 0 0 1

表3.2 第一次旋转后坐标

X Y Z 0 0 0 0.9888 0 -0.1494 0 1.0000 0 0.1494 0 0.9888 X Y Z 0 0 0 表3.3 第二次旋转后坐标 0.0699 0.9975 -0.9863 0.0707 -0.1494 0 0.0106 -0.1491 0.9888 《数值计算方法》实验报告 12

表3.4 第三次旋转后坐标

X Y Z 0 0 0 0.0699 0.9555 -0.2864 0.9975 -0.0640 0.0302 0.0106 -0.2878 -0.9576

P108.1

算法：

（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）.

流程图：

《数值计算方法》实验报告

output err=abs(norm(X'-P)); relerr=err/(norm(X)+eps);P=X'; X(j)=(B(j)-A(j,j-1)*X(j-1)'-A(j,j+1)*P(j+1))/A(j,j) X(N)=(B(N)-A(N,N-1)*(X(N-1))')/A(N,N) X(1)=(B(1)-A(1,2)*P(2))/A(1,1) j=1 k=k+1 k>max1 Y N=length(B); k=1 start InputA,B,P,delta,max1 N j=j+1 j>N Y N Y j==1 Y j==N N N N err

P109.7 算法：

（1）输入A,B，令[N N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B]; 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))，输出结果。

《数值计算方法》实验报告

end output X=backsub(Aug(1:N,1:N),Aug(1:N,N+1)) m=Aug(k,p)/Aug(p,p); N Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1) k=p+1 Aug(p,p)==0 p=p+1 start [N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B]; p=1 p>N-1 Y N [Y,j]=max(abs(Aug(p:N,p))); C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C; Y N Y k>N k=k+1

《数值计算方法》实验报告

P130.4

function X=gseid1(A,B,P,delta,max1) %input--A is an N*N nonsingular matrix % -B is an N*1 matrix % -P is an N*1 matrix

% -delta is the tolerance for P

% -max1 is the maximum number of iterations

%output-X is an N*1 matrix:the gauss-seidel approximation to the solution

%of AX=B

N=length(B);

for k=1:max1 for j=1:N if j==1

X(1)=(B(1)-A(1,2:3)*P(2:3))/A(1,1); elseif j==2

X(2)=(B(2)-A(2,1)*X(1)'-A(2,3:4)*P(3:4))/A(2,2); elseif 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);