数值分析-列主元高斯消去法
- 格式:docx
- 大小:190.50 KB
- 文档页数:1
一.实验要求二.算法描述三.实验代码//列主元高斯消元法#include<iostream>#include<cmath>using namespace std;double A[100][100];//系数矩阵double b[100];//右端项double x[100];//结果double e;//用于溢出控制int n;//矩阵n*n大小int main(){cout << "----列主元高斯消元法----\n";cout << "请输入矩阵的大小n及e:"; cin >> n >> e;//1.输入A,b,ecout << "请输入系数矩阵:" << endl;for (int i = 1; i <= n; i++)for (int j = 1; j <=n; j++)cin >> A[i][j];cout << "请输入右端矩阵:" << endl;for (int i = 1; i <= n; i++)cin >> b[i];for (int k = 1; k <= n - 1; k++)//2.选主元及消元{double T = A[k][k]; int ik = k;for(int l = k + 1; l <= n; l++)if(A[l][k]>T){ T = A[l][k]; ik = l; }if (T < e){ cout << "求解失败." << endl; exit(0);}if (ik != k){double temp;for (int l = 1; l <= n; l++)//交换ik行和k行{temp = A[ik][l]; A[ik][l] = A[k][l]; A[k][l] = temp;}temp = b[ik]; b[ik] = b[k]; b[k] = temp;//交换b_ik 和bk}for (int i = k + 1; i <= n; i++){T = A[i][k] / A[k][k];b[i] = b[i] - T*b[k];for (int j = k + 1; j <= n; j++){A[i][j] = A[i][j] - T*A[k][j];}}}//3.回代if (abs(A[n][n]) <= e){ cout << "求解失败." << endl; exit(0); }x[n] = b[n] / A[n][n];for (int i = n - 1; i >= 1; i--){double sum = 0;for (int j = i + 1; j <= n; j++)sum += A[i][j] * x[j];x[i] = (b[i] - sum) / A[i][i];}cout << "结果:" << endl;//打印xifor (int i = 1; i <= n; i++)cout << x[i] << " ";cout << endl;return 0;}//直接LU分解法#include<iostream>#include<cmath>using namespace std;double A[100][100];//系数矩阵double b[100];//右端项double u[100][100];//u矩阵double l[100][100];//l矩阵double y[100];//中间变量double x[100];//结果double e;//用于溢出控制int n;//矩阵n*n大小int main(){cout << "----直接LU分解法----\n";cout << "请输入矩阵的大小n及e:"; cin >> n >> e;//1.输入A,b,ecout << "请输入系数矩阵:" << endl;for (int i = 1; i <= n; i++)for (int j = 1; j <= n; j++)cin >> A[i][j];cout << "请输入右端矩阵:" << endl;for (int i = 1; i <= n; i++)cin >> b[i];//2.计算u_kj,l_ikfor (int k = 1; k <= n; k++){for (int j = k; j <= n; j++){double sum = 0;for(int m = 1; m <= k - 1; m++)sum += l[k][m] * u[m][j];u[k][j] = A[k][j] - sum;if (abs(u[k][k])<e){ cout << "求解失败." << endl;exit(0); }for (int i = k + 1; i <= n; i++){sum = 0;for (int m = 1; m <= k - 1; m++)sum += l[i][m] * u[m][k];l[i][k] = (A[i][k] - sum) / u[k][k];}}//3.求解LY=by[1] = b[1];for (int i = 2; i <= k; i++){double sum = 0;for(int j = 1; j <= i - 1; j++)sum += l[i][j] * y[j];y[i] = b[i] - sum;}//4.求解UX=Yx[n] = y[n] / u[n][n];for (int i = n - 1; i >= 1; i--){double sum = 0;for(int j = i + 1; j <= n; j++)sum += u[i][j] * x[j];x[i] = (y[i] - sum) / u[i][i];}}cout << "结果:" << endl;//5.打印xifor (int i = 1; i <= n; i++)cout << x[i] << " ";cout << endl;return 0;}四.实验结果。
课程设计任务书前 言回顾普通解方程组的方法,一般都是先逐个削去未知变量,最终得到只有一个未知变量的方程,解之,把得到的值回代到消去变量过程中得到的方程组,逐个求出未知变量。
这种解线性方程组的基本方法就是这里要介绍的高斯消去法。
数学上,高斯消元法(或译:高斯消去法),是线性代数中的一个算法,可用来为线性方程组求解,求出矩阵的秩,以及求出可逆方阵的逆矩阵。
当用于一个矩阵时,高斯消元法会产生出一个“行梯阵式”。
高斯消元法可以用在电脑中来解决数千条等式及未知数。
高斯消元法可以用来找出一个可逆矩阵的逆矩阵。
用关联矩阵表述网络拓扑结构,并根据厂站拓扑结构和网络拓扑结构等概念简化了电力系统的拓扑结构。
根据广义乘法和广义加法的运算规则,将改进的高斯消元算法应用于电力系统拓扑结构分析中,并引入稀疏、分块处理等技术提高了上述拓扑分析的效率。
采用上述高斯消元算法对山东电网220kV 以上的变电站进行拓扑结构分析,结果表明了运用该高斯消元法进行网络拓扑分析的正确性和有效性。
用列主元素法,选取每列的绝对值最大的元素作为消去对象并作为主元素。
然后换行使之变到主元位子上,在进行消元计算。
设)()(k k b X A ,确定第k 列主元所在位置k i ,在交换k i 行和k 行后,在进行消元,并用MATLAB 软件进行求解。
目录摘要......................................................................................... 错误!未定义书签。
第1章绪论 ........................................................................... 错误!未定义书签。
第2章高斯消元法的算法描述 (2)2.1高斯消元法的原理概述 (2)2.1.1高斯消元法的消元过程 (2)2.1.2高斯消元法的回带过程 (3)2.1.3高斯消元法的复杂度分析 (4)2.2列主高斯消元法原理简介 (5)2.2.1列主高斯消元法的消元过程 (6)2.2.2列主高斯消元法的回带过程 (6)2.2.3列主高斯消元法的算法描述 (6)第3章高斯消元法的物理应用 (9)3.1电网模型的描述 (9)3.2电网模型的问题分析 (9)3.3求解计算 (11)参考文献 (13)摘 要用列主元素高斯消去法法,选取每列的绝对值最大的元素作为消去对象并作为主元素。
.精品精品第5章复习与思考题1、用高斯消去法为什么要选主元?哪些方程组可以不选主元?答:使用高斯消去法时,在消元过程中可能出现0k kka = 的情况,这时消去法无法进行;即时主元素0kkk a ≠,但相对很小时,用其做除数,会导致其它元素数量级的严重增长和舍入误差的扩散,最后也使得计算不准确。
最后也使得计算不准确。
因此高斯消去法需要选主元,因此高斯消去法需要选主元,因此高斯消去法需要选主元,以保证计算的进行和计以保证计算的进行和计算的准确性。
当主对角元素明显占优(远大于同行或同列的元素)时,当主对角元素明显占优(远大于同行或同列的元素)时,可以不用选择主元。
可以不用选择主元。
可以不用选择主元。
计算时一般选计算时一般选择列主元消去法。
2、高斯消去法与LU 分解有什么关系?用它们解线性方程组Ax = b 有何不同?A 要满足什么条件?答:高斯消去法实质上产生了一个将A 分解为两个三角形矩阵相乘的因式分解,其中一个为上三角矩阵U ,一个为下三角矩阵L 。
用LU 分解解线性方程组可以简化计算,减少计算量,提高计算精度。
A 需要满足的条件是,顺序主子式(需要满足的条件是,顺序主子式(1,21,21,2,…,,…,,…,n-1n-1n-1)不为零。
)不为零。
3、楚列斯基分解与LU 分解相比,有什么优点?楚列斯基分解是LU 分解的一种,当限定下三角矩阵L 的对角元素为正时,的对角元素为正时,楚列斯基分解具楚列斯基分解具有唯一解。
4、哪种线性方程组可用平方根法求解?为什么说平方根法计算稳定?具有对称正定系数矩阵的线性方程可以使用平方根法求解。
平方根法在分解过程中元素的数量级不会增长,平方根法在分解过程中元素的数量级不会增长,切对角元素恒为正数,切对角元素恒为正数,因此,是一个稳定的算法。
5、什么样的线性方程组可用追赶法求解并能保证计算稳定? 对角占优的三对角方程组6、何谓向量范数?给出三种常用的向量范数。
一、 列主元素Gauss 消去法、Jacobi 迭代法原理及计算方法1. 列主元素Gauss 消去法:1.1 Gauss 消去法基本原理设有方程组Ax b =,设A 是可逆矩阵。
高斯消去法的基本思想就是将矩阵的初等行变换作用于方程组的增广矩阵[]B A b = ,将其中的A 变换成一个上三角矩阵,然后求解这个三角形方程组。
1.2 列主元Gauss 消去法计算步骤将方程组用增广矩阵[]()(1)ijn n B A b a ⨯+== 表示。
1). 消元过程对1,2,,1k n =-(1) 选主元,找{},1,,k i k k n ∈+ 使得 ,max k i k ik k i na a ≤≤= (2) 如果,0k i k a =,则矩阵A 奇异,程序结束;否则执行(3)。
(3) 如果k i k ≠,则交换第k 行与第k i 行对应元素位置,k kj i j a a ↔,,,1j k n =+ 。
(4) 消元,对,,i k n = ,计算/,ik ik kk l a a =对1,,1j k n =++ ,计算.ij ij ik kj a a l a =-2). 回代过程(1) 若0,nn a =则矩阵奇异,程序结束;否则执行(2)。
(2) ,1/;n n n nn x a a +=对1,,2,1i n =- ,计算,11/n i i n ij j ii j i x a a x a +=+⎛⎫=- ⎪⎝⎭∑2. Jacobi 迭代法2.1 Jacobi 迭代法基本原理Jacobi 迭代法的基本思想是对n 元线性方程组b Ax =,.,n n R b R A ∈∈将其变形为等价方程组f Bx x +=,其中.,,n n n n R x R f R B ∈∈∈⨯B 成为迭代矩阵。
从某一取定的初始向量)0(x 出发,按照一个适当的迭代公式 ,逐次计算出向量f Bx x k k +=+)()1( ( 1,0=k ),使得向量序列}{)(k x 收敛于方程组的精确解.(1)输入1,,,,)0(=k n xb A ε,. (2) )(1,1)0()1(∑≠=-=n j i i j ij i iii x a b a x )1,0(n i = (3)判断 ε≤--≤≤)0()1(10max i i n i x x ,若是,输出1)1(2)1(1,,n x x x ,若否,置1+=k k ,)1()0(i i x x =,)2,1(n i =。
数值分析1顺序消去法、列主元、列主元Gauss-Jordan消去法function x = Gauss (A, b)% 求解方程组的Gauss消去法,调用方法为% x = Gauss (A, b)% 其中% A 为方程组的系数矩阵,b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);if n~=merror ('% 系数矩阵必须为方的!');endif m~=nberror ('% b 的维数与方程组的行数不匹配!'); endfor k = 1:n-1% 消元过程for i = k+1:nm = A (i,k)/A(k,k);for j = k+1:nA (i,j) = A (i,j)-m*A (k,j);endb (i) = b (i)-m*b (k);endendx=zeros (size (b));for k = n:-1:1for j = k+1:nb (k) = b (k)-A (k,j)*x (j);endx (k) = b (k)/A(k,k);endendfunction x = Gauss_Elim (A, b)% 求解方程组的列主元Gauss消去法,调用方法为% x = Gauss_Elim (A, b)% 其中% A 为方程组的系数矩阵,b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);error ('% 系数矩阵必须为方的!');endif m~=nberror ('% b 的维数与方程组的行数不匹配!');endfor k = 1:n-1% 选主元a_max = 0;for i = k:nif abs (A (i,k))>a_maxa_max = A (i,k); r=i;endendif abs(a_max)<1e-15error ('% 系数矩阵奇异,无法求解方程组!');end% 交换两行if r>kfor j = k:nz=A (k,j); A (k,j)=A (r,j);A (r,j)=z;endz=b (k);b (k)=b (r);b (r)=z;end% 消元过程for i = k+1:nm = A (i,k)/A(k,k);for j = k+1:nA (i,j) = A (i,j)-m*A (k,j);endb (i) = b (i)-m*b (k);endend% 回代过程if abs (A (n,n))<1e-15error ('% 系数矩阵奇异,无法求解方程组!'); endx=zeros (size (b));for k = n:-1:1for j = k+1:nb (k) = b (k)-A (k,j)*x (j);endx (k) = b (k)/A(k,k);endendfunction x = Gauss_Jordan (A, b)% 求解方程组的列主元Gauss-Jordan消去法,调用方法为% x = Gauss_Jordan (A, b)% 其中% A 为方程组的系数矩阵,b为方程组的右端项% x 为方程组的解[n,m] = size (A); nb = length (b);error ('% 系数矩阵必须为方的!');endif m~=nberror ('% b 的维数与方程组的行数不匹配!'); endfor k = 1:n% 选主元a_max = 0;for i = k:nif abs (A (i,k))>a_maxa_max = A (i,k); r=i;endendif abs(a_max)<1e-15error ('% 系数矩阵奇异,无法求解方程组!'); end% 交换两行if r>kfor j = k:nz=A (k,j); A (k,j)=A (r,j);A (r,j)=z;endz=b (k);b (k)=b (r);b (r)=z;end% 消元计算b (k) = b (k)/A(k,k);for j = k+1:nA (k,j) = A (k,j)/A(k,k);endfor i=1:nfor j=k+1:nA (i,j) = A (i,j)-A (i,k)*A (k,j); endb (i)=b (i)-A (i,k)*b (k); endendendx = b; % 输出bend。
实验四:列组元消去法一、目的1)熟悉列主元高斯消元法解线性方程组的算法2)掌握列主元高斯消去法的编程二、实验原理列主元素消去法是为控制舍入误差而提出来的一种算法,在Gauss消去法的消元过程中,若出现a=0,则消元无法进行,即使其不为0,但很小,把它作为除数,就会导致其他元素量级的巨大增长和舍入误差的扩散,最后使计算结果不可靠.使用列主元素消去法计算,基本上能控制舍入误差的影响,并且选主元素比较方便.三、运行结果四、代码using System;using System.Collections.Generic;using System.Linq;using System.Text;namespace高斯{class Program{static double[] Gause(double[,] a, int n){int i, j, k;int rank, columm;double temp, l, s, mx;double[] x = new double[n];for (i = 0; i <= n - 2; i++){mx = Math.Abs(a[i, i]);rank = i;columm = i;for (j = i + 1; j <= n - 1; j++) //选主元if (Math.Abs(a[j, i]) > mx){mx = Math.Abs(a[j, i]);rank = j;columm = i;}for (k = 0; k <= n; k++) //主元行变换{temp = a[i, k];a[i, k] = a[rank, k];a[rank, k] = temp;} //消元for (j = i + 1; j <= n - 1; j++){l = a[j, i] / a[i, i];for (k = i; k <= n; k++)a[j, k] = a[j, k] - l * a[i, k];}}x[n - 1] = a[n - 1, n] / a[n - 1, n - 1]; //回代方程求解x for (i = n - 2; i >= 0; i--){s = 0;for (j = i + 1; j <= n - 1; j++)s = s + a[i, j] * x[j];x[i] = (a[i, n] - s) / a[i, i];}return x;}static void Main(string[] args){double[,] a = new double[4, 5] { { 10, -7, 0, 1, 8 }, { -3, 2.099999, 6, 2, 5.900001 }, { 5, -1, 5, -1, 5 }, { 2, 1, 0, 2, 1 } };int n = 4;double[] x = new double[n];x = Gause(a, n);Console.WriteLine("高斯消去法方程:");for (int i = 0; i < n; i++){for (int j = 0; j < n; j++)Console.Write(a[i, j].ToString() + " ");Console.WriteLine();}Console.WriteLine("线性方程组的解:");for (int i = 0; i <= n - 1; i++)Console.Write("x" + (i + 1).ToString() + "=" +x[i].ToString() + " ");Console.WriteLine();Console.ReadLine();}}}四、分析通过本次实验的学习,学会根据算法编写基本的相关程序,虽然此次程序模板由老师给予,但认真阅读理解程序有助于今后的学习,再利用计算机中的C语言对高斯列主元消去法可以快速得到线性方程组的解,由简单的线性方程组可以推广到一般n阶线性方程组,这对如何利用高斯列主元消去法解决实际问题有了一定的经验。
§2高斯主元素消去法⎪⎩⎪⎨⎧=++-=++=++00.357.404.100.200.224.563.200.100.100.200.10120.0321321321x x x x x x x x x 解:clear alla=[0.0120 1.00 2.00;1.00 2.63 5.24;-2.00 1.04 4.57]; b=[1.00;2.00;3.00];x=a\b方程组的三位有效数字的解:Tx )266.0,476.0,645.0(*-=Gauss 消去法求解(取三位有效数字):[]⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡------−−−→−⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-==-==00.300.5003.811627.80000.100.200.10120.016432916603.811627.80000.100.200.10120.000.357.404.100.200.224.563.200.100.100.200.10120.006.21673.83323121l l l b A 解出Tx )60.0,197.0,25.0(--≈。
【注】1)设Ax=b,其中A 为n 阶非奇异矩阵,可以应用高斯消元法。
2)消元过程中,即使0)(≠k kk a ,用其作除数)/()()(k kk k ik ik a a l =会导致计算中间结果数量级严重增长和舍入误差的累积、扩大,最后使得计算结果不可靠。
3)应避免采用绝对值很小的主元素)(k kk a ;对一般的系数矩阵,最好保持乘数1≤ik l ,因此,在高斯消去法中应引进选主元技巧,以便减少计算过程中舍入误差对求解的影响。
clear alla=[0.0120 1.00 2.00;1.00 2.63 5.24;-2.00 1.04 4.57]; b=[1.00;2.00;3.00];x_value=vpa(a\b,15)%10位有效数字的近似解a=[a,b];eps=1e-6;[n,m]=size(a);Gauss,x=vpa(x,15) %对比高斯消去法的结果一、列主元素消去法基本思想:在每轮消元之前,选列主元素(绝对值最大的元素),使乘数(即消元因子)1≤ik l步骤:设已进行k-1轮消元,得矩阵⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=)()()()()2(2)2(2)2(22)1(1)1(1)1(12)1(11)(k nn k nkk kn k kkn kn kk a a a a a a a a a a a A一落千丈 1 23S1:选列主元素: )()(0max k ik ni k k k i a a ≤≤= (1)S2:换行:如果)(0k k i a →0,则方程组解不唯一,停止运算; 否则,如果i0=k , 则可进行下一轮消元;如果k i ≠0,则r i0 r k ,然后进行下一轮消元。
数值分析实验报告(1)学院:信息学院班级:计算机0903班姓名:***学号:********课题一A.问题提出给定下列几个不同类型的线性方程组,请用适当的方法求解线性方程组1、设线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------------------1368243810041202913726422123417911101610352431205362177586832337616244911315120130123122400105635680000121324⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x =⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2119381346323125 x *= ( 1, -1, 0, 1, 2, 0, 3, 1, -1, 2 )T2、设对称正定阵系数阵线方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----------------------19243360021411035204111443343104221812334161206538114140231212200420424⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡87654321x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡---4515229232060 x * = ( 1, -1, 0, 2, 1, -1, 0, 2 )T3、三对角形线性方程组⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡------------------4100000000141000000001410000000014100000000141000000001410000000014100000000141000000001410000000014⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡10987654321x x x x x x x x x x = ⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡----5541412621357 x *= ( 2, 1, -3, 0, 1, -2, 3, 0, 1, -1 )TB.(1)对上述三个方程组分别用Gauss 顺序消去法与Gauss 列主元消去法;平方根 与改进平方根法;追赶法求解(选择其一) (2)编写算法通用程序(3)在应用Gauss 消去时,尽可能利用相应程序输出系数矩阵的三角分解式C.(1)通过该课题的程序编制,掌握模块化结构程序设计方法 (2)掌握求解各类线性方程组的直接方法,了解各种方法的特点 (3)体会高斯消去法选主元的必要性 实验步骤:(高斯消去法,列主元,LU )1顺序高斯消去法2.LU 分解法3.列主元高斯消去法(如下图)(1)高斯消去法运行结果如下(2)对方程的系数矩阵进行LU分解并求出方程组的解(3)列主元高斯消去法实验体会总结:利用gauss消去法解线性方程组的时候,如果没有经过选主元,可能会出现数值不稳定的现象,使得方程组的解偏离精确解。
《计算方法》实验指导书 实验三、高斯消元法和列主消元法一、实验目的:1. 通过matlab 编程解决高斯消元发和列主消元发来解方程组的问题, 加强编程能力和编程技巧,要熟练应用matlab 程序来解题,练习从数值分析的角度看问题进而来解决问题。
更深一步体会这门课的重要性,练习动手能力,同时要加深对数值问题的理解,要熟悉matlab 编程环境。
二、实验要求:用matlab 编写代码并运行高斯消元法和列主消元发来解下面的方程组的问题,并算出结果。
三、实验内容:用高斯消元法和列主消元法来解题。
1.实验题目:用高斯消元法和列主消元法来解下列线性方程组。
⎪⎪⎩⎪⎪⎨⎧−=+−−−=+−−=+−−=−+−.142,16422,0,13143214321432432x x x x x x x x x x x x x x x 2.实验原理高斯消元法:就是把方程组变成上三角型或下三角形的解法。
上三角形是从下往上求解,下三角形是从上向下求解,进而求得结果。
而列主消元法是和高斯消元法相类似,只不过是在开始的时候找出x1的系数的最大值放在方程组的第一行,再化三角形再求解。
3.设计思想高斯消元法:先把方程组的第一行保留,再利用第一行的方程将其余几行的含有x1的项都消去,再保留第二行,同理利用第二行的方程把第二行以下的几行的含有x2项的都消去,以此类推。
直到最后一行只含有一个未知数,化为上三角形,求得最后一行的这个未知数的值,再回带到倒数第二个方程求出另一个解,再依次往上回带即可求出这个方程组的值。
而列主消元法与高斯消元法类似,只不过在最开始时找出x1项系数的最大值与第一行交换再进行与高斯算法相似的运算来求出方程组的解。
4.源代码高斯消元法的程序:f unction [RA,RB,n,X]=gaus(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endend在工作窗口输入程序:A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1];b=[1;0; -1;-1]; [RA,RB,n,X] =gaus (A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.运行结果为:RA =4RB =4n =4X =-0.50000.5000.列主消元发的程序:function [RA,RB,n,X]=liezhu(A,b)B=[A b]; n=length(b); RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意:因为RA~=RB,所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意:因为RA=RB=n,所以此方程组有唯一解.')X=zeros(n,1); C=zeros(1,n+1);for p= 1:n-1[Y,j]=max(abs(B(p:n,p))); C=B(p,:);B(p,:)= B(j+p-1,:); B(j+p-1,:)=C;for k=p+1:nm= B(k,p)/ B(p,p);B(k,p:n+1)= B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意:因为RA=RB<n,所以此方程组有无穷多解.')endend在工作窗口输入程序:A=[1 -1 1 -3; 0 -1 -1 1;2 -2 -4 6;1 -2 -4 1];b=[1;0; -1;-1]; [RA,RB,n,X]=liezhu(A,b)请注意:因为RA=RB=n,所以此方程组有唯一解.运行结果为:RA =4RB =4n =4X =-0.50000.5000实验体会:通过这次实验我了解了高斯消元法和列主消元方法的基本思想,虽然这两个程序的编写是有点困难的,但运行起来还是比较容易的,解决了不少实际问题的计算。
课题名称:课题一解线性方程组的直接方法解决的问题:给定三个不同类型的线性方程组,用适当的直接法求解。
采用的数值方法:对第一个普通的线性方程组,采用了高斯顺序消去法和高斯列主元消去法。
对第二个正定线性方程组,采用了平方根法。
对第三个三对角线性方程组,采用了追赶法。
算法程序:(1)普通的线性方程组①顺序消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0},{8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;int size=10;for(k=0; k<size-1; k++){if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k];for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j];}x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}②列主元消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0}, {8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;float max;int col;int size=10;for(k=0; k<size-1; k++){max=fabs(A[k][k]);col=k;for(i=k; i<size; i++){if(max<fabs(A[i][k])) {max=fabs(A[i][k]); col=i;}}for(j=k; j<size; j++){temp=A[col][j];A[col][j]=A[k][j];A[k][j]=temp;}temp=b[col];b[col]=b[k];b[k]=temp;if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k]; for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j]; }x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}(2)对称正定线性方程组平方根法:#include <stdio.h>#include <math.h>#define n 8int main(void){float A[8][8]={{4,2,-4,0,2,4,0,0},{2,2,-1,-2,1,3,2,0},{-4,-1,14,1,-8,-3,5,6},{0,-2,1,6,-1,-4,-3,3},{2,1,-8,-1,22,4,-10,-3},{4,3,-3,-4,4,11,1,-4},{0,2,5,-3,-10,1,14,2},{0,0,6,3,-3,-4,2,19}};float g[8][8]= {0};float b[8]= {0,-6,6,23,11,-22,-15,45}; float x[8]= {0};float y[8]= {0};int k,m,i,sq;for(k=0; k<n; k++){float p=0,q=0,s=0;for(m=0; m<=k-1; m++){p=p+A[k][m]*A[k][m];}g[k][k]=sqrt(A[k][k]-p);A[k][k]=g[k][k];for(i=k+1; i<n; i++){q=0;for(m=0; m<=k-1; m++){q=q+A[i][m]*A[k][m];}g[i][k]=(A[i][k]-q)/A[k][k]; A[i][k]=g[i][k];}s=0;for(m=0; m<=k-1; m++){s=s+A[k][m]*y[m];}y[k]=(b[k]-s)/A[k][k];}x[n-1]=y[n-1]/A[n-1][n-1];for(k=n-2; k>=0; k--){float sum=0;for(m=k+1; m<n; m++){sum=sum+A[m][k]*x[m];}x[k]=(y[k]-sum)/A[k][k];}for(sq=0; sq<n; sq++){printf("%f ",x[sq]);}return 0;}(3)三对角线性方程组追赶法#include <stdio.h>#include <math.h>#define n 10int main(void){float a[10]={4,4,4,4,4,4,4,4,4,4};float c[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1};float d[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1}; float b[10]={7,5,-13,2,6,-12,14,-4,5,-5}; float x[10]={0};float y[10]={0};float arf[10]={0};float bt[9]={0};arf[0]=a[0];int i;for(i=0;i<n-1;i++){bt[i]=c[i]/arf[i];arf[i+1]=a[i+1]-d[i+1]*bt[i];//printf("%f %f \n",bt[i],arf[i+1]); }y[0]=b[0]/arf[0];//printf("%f\n",y[0]);for(i=1;i<n;i++){y[i]=(b[i]-d[i]*y[i-1])/arf[i];}//printf("%f\n",y[1]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--){x[i]=y[i]-bt[i]*x[i+1]; }for(i=0;i<n;i++)printf("%lf ",x[i]);return 0;}数值结果:(1) 普通的线性方程组①顺序消去法②列主元消去法(2) 对称正定线性方程组平方根法:(3)三对角线性方程组追赶法:对实验计算结果的讨论和分析:(1) 普通的线性方程组①顺序消去法x1~x10的绝对误差:0.000001,-0.000001,0.000001,0,0.000001,0,0.000002,0,0,0x1~x10的相对误差:0.000001,0.000001,-1,0,0.0000005,0,0.00000067,0,0,0误差很小,基本可以忽略。