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

数值计算方法 实验报告1 线性方程组线性方程组 AX BAX 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 2 实验内容实验内容 P93 1 单位立方体位于第一卦限 一个顶点在原点 首先 以角度 6 沿 y 轴旋转立方体 然后 再以角度 4 沿 z 轴旋转立方体 求旋转后立方体的 8 个顶点的坐标 并与例 3 10 的结果比 较 它们的区别是什么 试通过矩阵一般不满足交换律的事实进行解释 使用 plot3 命令画 出 3 个图形 P93 2 设单位立方体位于第一卦限 其中一个顶点位于坐标原点 首先以角度 12 沿 x 轴旋转 立方体 然后再以角度 6 沿 z 轴旋转立方体 求旋转后立方体的 8 个顶点的坐标 使用 plot3 画出这 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 许多科学应用包含的矩阵带有很多零 在实际情况中很重要的三角形线性方程组有如下 形式 1 1121 d xc xb 数值计算方法 实验报告 2 1 122232 a xd xc xb 2233343 a xd xc xb 221111NNNNNNN axdxcxb 11NNNNN axd xb 构造一个程序求解三角形线性方程组 可假定不需要变换 而且可用第 k 行消去第 k 1 行的 k x p109 7 下面的习题虽然是针对 3x3 维矩阵的 但其概念可用于 NxN 维矩阵 如果矩阵 A 非奇 异 则 1 A 存在 而且 1 AAI 设 1 C 2 C 3 C是 1 A 的列 而 1 E 2 E 3 E是 I 的列 方程 1 AAI 可表示为 123123 A C C CE E E 则上式等价于三个线性方程组 11 ACE 22 ACE 33 ACE 这样求 1 A 等价于求解三个线性方程组 使用程序 2 2 或上题中的程序求解下面每个矩阵的逆 通过计算 1 AA 和使用命令 inv A 检 查答案 并解释可能的差异 a 201 325 110 b 16120240140 120120027001680 240270064804200 140168042002800 120 3 修改程序 3 3 使得它可以通过重复求解 N 个线性方程组 JJ ACE 其中 J 1 2 N 来得到 1 A 则 1212 NN A CCCEEE 而且 1 12 N ACCC 保证对 LU 分解只计算一次 p120 4 基尔霍夫电压定律说明电路网络中任意单向闭路的电压和为零 现有如下电路线性方程 数值计算方法 实验报告3 134132431 312352532 41524563 0 RRRIR IR IE R IRRRIR IE R IR IRRRI 如果 12345612 12345612 12345612 1 1 2 1 2 4 23 29 1 0 75 1 21 4 12 21 5 1 2 4 3 1 5 41 38 a RRRRRREE b RRRRRREE c RRRRRREE 使用程序求解电流 123 I II p129 4 利用高斯 赛德尔迭代法求解下列带状方程 123 1234 12345 23456 4647484950 47484950 484950 1225 21225 21225 21225 21225 21225 2125 xxx xxxx xxxxx xxxxx xxxxx xxxx xxx 3 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 100 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 100 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 plo t3 W d 1 W d 2 W d 3 draw now 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 Y Y N N i i 1 i i 1 数值计算方法 实验报告5 0 0 5 1 0 0 5 1 0 0 5 1 1 0 1 0 0 5 1 0 1 2 2 0 2 2 0 2 0 1 2 x y z 图1 1 左上 初始立方体 图1 2 右上 沿y轴旋转 6 y VRU 图1 3 左下 沿z轴旋转 4 x WRV 表1 1 第一次旋转后立方体的坐标 表 1 2 第二次旋转后立方体的坐标 X00 86600 866000 50000 50001 36601 3660 Y001 00001 000 0 1 0000001 0000 Z0 0 5000 0 500000 86600 86600 36600 3660 X00 6124 0 0947 0 7071 0 3536 0 35360 96590 2588 Y00 61241 31950 70711 06070 35360 96591 6730 Z0 0 5000 0 5000 00 86600 86600 36600 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 数值计算方法 实验报告7 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 100 r1 1 0 0 0 cos i pi 1200 sin i pi 1200 0 sin i pi 1200 c os i pi 1200 U X r1 plot3 U d 1 U d 2 U d 3 drawnow i 0 i 100 r2 cos i pi 600 sin i pi 600 0 sin i pi 600 cos i pi 6 00 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 Y Y N N i i 1 i i 1 数值计算方法 实验报告 8 0 0 5 1 0 0 5 1 0 0 5 1 0 0 5 1 1 0 1 0 1 2 1 0 1 2 0 2 0 1 2 x y z 图2 1 左上 单位立方体 图2 2 右上 第一次旋转 图2 3 左下 第2次旋转 表2 1 第一次旋转后立方体坐标 表2 2 第二次旋转后立方体坐标 X01001101 Y000 9659 0 25880 9659 0 25880 70710 7071 Z000 25880 96590 25880 96591 22471 2247 X00 8660 0 48300 12940 38310 9954 0 35360 5125 Y00 50000 8365 0 2241 1 33650 27590 61241 1124 Z000 25880 96590 25880 96591 22471 2247 数值计算方法 实验报告9 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 xlabel x ylabel y zlabel z view 3 rotate3d 数值计算方法 实验报告 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 100 r1 1 0 0 0 cos i pi 1200 sin i pi 1200 0 sin i pi 1200 c os i pi 1200 U X r1 plot3 U d 1 U d 2 U d 3 drawnow i 0 i 100 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 Y Y N N i i 1 i i 1 i 0 i 100 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 N N i i 1 数值计算方法 实验报告11 0 0 5 1 0 0 5 1 0 0 5 1 1 0 1 0 0 5 1 1 0 1 0 0 1 0 2 2 0 2 1 0 1 0 0 1 0 2 1 0 1 1 0 1 x y z 图1 1 左上 正四面体 图1 2 右上 第一次旋转 图3 3 左下 第二次旋转 图3 4 右下 第4次旋转 表3 1 四面体坐标 表3 2 第一次旋转后坐标 表3 3 第二次旋转后坐标 X0100 Y0010 Z0001 X00 988800 1494 Y001 00000 Z0 0 149400 9888 X00 06990 99750 0106 Y0 0 98630 0707 0 1491 Z0 0 149400 9888 数值计算方法 实验报告 12 表3 4 第三次旋转后坐标 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 delta relerrmax1 j 1 j N j 1 j N err delta rel errN 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 数值计算方法 实验报告15 start N N size A X zeros N 1 C zeros 1 N 1 A ug A B p 1 Y j max abs Aug p N p C Aug p Aug p Aug j p 1 Aug j p 1 C p N 1 Aug p p 0 k p 1 k N m Aug k p Aug p p Aug k p N 1 Aug k p N 1 m Aug p p N 1 output end X backsub Aug 1 N 1 N Aug 1 N N 1 Y N Y N Y N k k 1 p p 1 数值计算方法 实验报告 16 a b 1 00 20 4 1 00 21 4 1 00 40 8 1 00000 50000 33330 2500 0 50000 33330 25000 2000 0 33330 25000 20000 1667 0 25000 20000 16670 1429 P120 3 修改程序3 3 主要变动为列矩阵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 输出结果 数值计算方法 实验报告17 流程图 start inputA E p 1 N 1 N N size A X zeros N 1 Y zeros N 1 C zeros 1 N t 1 N p 1 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 A p p 0 k p 1 k N 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 J 1 J N k 2 Y 1 E t 1 J k N Y k E t k J A k 1 k 1 Y 1 k 1 X N Y N A N N k N 1 k 1 Y k k 1 N Y N Y N k k 1 YJ J 1 Yk k 1 N Yk k 1 N N 数值计算方法 实验报告 18 P120 4 算法 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 输出结果 X k Y k A k k 1 N X k 1 N A k k m 1 N J X 1 N output end 数值计算方法 实验报告19 流程图 数值计算方法 实验报告 20 start inputA E N N size A X zeros N 1 Y zeros N 1 C zeros 1 N t 1 N p 1 p 1 N 1 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 A p p 0 k p 1 k N 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 2 Y 1 E t 1 J k N Y k E t k J A k 1 k 1 Y 1 k 1 X N Y N A N N k N 1 output end k 1 X k Y k A k k 1 N X k 1 N A k k k k 1Y N Y N Yk k 1 N Y k k 1 N Yk k 1 N 数值计算方法 实验报告21 表 3 电流 123 I II P130 4 算法 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 delta relerrmax1 j N j 1 j 1 X 1 B 1 A 1 2 P 2 A 1 1 j N X N B N A N N 1 X N 1 A N N err abs norm X P relerr err norm X eps P X err delta rel err delta output end j 2 X 2 B 2 A 2 1 X 1 A 2 3 4 P 3 4 A 2 2 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 X j B j A j j 1 X j 1 A j j 1 P j 1 A j j Y N Y N Y N Y N Y N Y N Y k k 1 j j 1 数值计算方法 实验报告23 表4 通过高斯赛德尔迭代得到的x的解 x 1 1 0 0 46 37955 23816 55 0 53 72846 05199 9656 0 50 90229 24601 3306 0 49 82216 34436 1741 0 49 89418 60239 7619 0 49 99853 51248 1308 0 50 00887 23890 1357 0 50 00153 18846 052 0 49 99947 93266 9753 0 49 99978 56913 4675 11 20 0 50 00001 08425 1992 0 50 00002 01576 6873 0 50 00000 22610 945 0 49 99999 86238 5722 0 49 99999 95873 979 0 50 00000 00530 2938 0 50 00000 00439 0388 0 50 00000 00023 9392 0 49 99999 99966 016 0 49 99999 99992 5573 21 30 0 50 00000 00001 7429 0 50 00000 00000 9169 0 50 49 99999 99999 9205 0 50 50 49 99999 99999 9212 0 50 50 00000 00000 9176 0 50 00000 00001 7445 31 40 0 49 99999 99992 555 0 49 99999 99966 017 0 50 00000 00023 939 0 50 00000 00439 0391 0 50 00000 00530 2936 0 49 99999 95873 979 0 49 99999 86238 5723 0 50 00000 22610 9451 0 50 00002 01576 6873 0 50 00001 08425 199 41 50 0 49 99978 56913 4675 0 49 99947 93266 9753 0 50 00153 18846 052 0 50 00887 23890 1357 0 49 99853 51248 1308 0 49 89418 60239 7619 0 49 82216 34436 1741 0 50 90229 24601 3306 0 53 72846 05199 9657 0 46 37955 23816 5501 4 4 结论结论 P93 1 2 3 1 利用矩阵乘的方法可对已知三维图形进行坐标变换 实现该图形在空间的旋转 2 矩阵乘一般不满足交换律 由于利用矩阵乘的方法对三维图形进行坐标变换 例如 V R1 U W R2 V 此时 故在进行多次旋转时 需要分步进行 不可一步2 1 WRRU 完成旋转 P108 1 7 通过对系数矩阵的增广矩阵进行高斯消元和回带容易得到线性方程组的解 同时 利用这 种方法可以求得矩阵的逆 P120 3 4 如果系数矩阵能分解为LU的形式 其中L为下三角矩阵 U为上三角矩阵 通过对系数矩 阵的分解再求解可应用简单的迭代进行求解x 数值计算方法 实验报告 24 P129 4 求解带状线性方程组的解可使用高斯 赛德尔迭代法 附件 代码 附件 代码 P93 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 for i 0 1 100 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 end for i 0 1 100 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 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 数值计算方法 实验报告25 p93 2 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 for i 0 1 100 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 end for i 0 1 100 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 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 p93 3 X zeros 4 3 数值计算方法 实验报告 26 X 2 7 12 1 d 1 2 3 4 1 3 2 4 for i 0 1 100 r1 cos i 0 15 100 sin i 0 15 100 0 0 1 0 sin i 0 15 100 0 cos i 0 15 100 U X r1 plot3 U d 1 U d 2 U d 3 drawnow end for i 0 1 100 r2 cos i 1 5 100 sin i 0 15 100 0 sin i 0 15 100 cos i 0 15 100 0 0 0 1 W U r2 plot3 W d 1 W d 2 W d 3 drawnow end for i 0 1 100 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 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 P108 1 function X gseid2 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 数值计算方法 实验报告27 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 P 2 A 1 1 elseif X N B N A N N 1 X N 1 A N N else X contains the kth approximations and P the k 1 st X j B j A j j 1 X j 1 A j j 1 P j 1 A j j end end err abs norm X P relerr err norm X eps P X if err delta relerr delta break end end p109 7 上三角变换 function X uptrbk A B input A is an N N nonsingular matrix B is an N 1 matrix output X is an N 1 matrix containing the solution to AX B initialize X and the temporary storage matrix C N N size A X zeros N 1 C zeros 1 N 1 Form the augmented matrix Aug A B Aug A B for p 1 N 1 Partial proviting for column p Y j max abs Aug p N p 数值计算方法 实验报告 28 interchange row p and j C Aug p Aug p Aug j p 1 Aug j p 1 C if Aug p p 0 A was singular No unique solution break end Elimination process for cloumn p for k p 1 N m Aug k p Aug p p Aug k p N 1 Aug k p N 1 m Aug p p N 1 end end back substitustion on U Y using program 3 1 X backsub Aug 1 N 1 N Aug 1 N N 1 回代 function X backsub A B input A is an n n upper triangular nonsinggular matrix B is an n 1 matrix output X is the solution to the linear system AX B Find the dimension of B and initialize X n length B 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 a 输入 A 2 0 1 3 2 5 1 1 0 B 1 0 0 X1 uptrbk A B B 0 1 0 X2 uptrbk A B B 0 0 1 X3 uptrbk A B X1 X2 X3 inv A b 输入 A 16 120 240 140 120 1200 2700 1680 240 2700 6480 4200 140 1680 4200 2800 B 1 0 0 0 数值计算方法 实验报告29 X1 uptrbk A B B 0 1 0 0 X2 uptrbk A B B 0 0 1 0 X3 uptrbk A B B 0 0 0 1 X4 uptrbk A B X1 X2 X3 X4 inv A p120 3 function m newlufact A E input A is an N N nonsingular matrix E is an N N eyes output m is an N N invmatrix of A initialize X Y the temporary storage matrix C and the row permutation information matrix t N N size A X zeros N 1 Y zeros N 1 C zeros 1 N t 1 N for p 1 N 1 find the privot row for column p max1 j max abs A p N p interchange row p and j 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 if A p p 0 A is singular No unique solution break end calculate multiplier and place in subdiagonal protion of A for k p 1 N 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 数值计算方法 实验报告 30 end end solve for Y for J 1 N Y 1 E t 1 J for k 2 N Y k E t k J A k 1 k 1 Y 1 k 1 end solve for X X N Y N A N N for k N 1 1 1 X k Y k A k k 1 N X k 1 N A k k end m 1 N J X 1 N end P120 4 function X lufact A B input A is an N N nonsingular matrix B is an N 1 matrix output X is an N 1 matrix containing the solution to AX B initialize X Y the temporary storage matrix C and the row permutation info