数值分析实验2_求解线性方程组直接法

《数值分析》课程教学大纲课程编号：07054352课程名称：数值分析英文名称：Numerical Analysis课程类型：学科基础课程要求：必修学时/学分：48/3 （讲课学时：40 上机学时：8）适用专业：计算机科学与技术；软件工程一、课程性质与任务“数值分析”是计算机科学与技术、软件工程等相关专业学生的学科基础课，也是其它理、工科专业本科生及研究生的必修或选修课。

数值分析是研究各种数学问题在计算机上通过数值运算，得到数值解答的方法和理论。

随着计算机系统能力的提高和新型数值软件的不断开发，无论在高科技领域还是在传统学科领域，数值分析的理论和方法的作用和影响巨大，是科学工作者和工程技术人员必备的基础知识和工具。

课程的任务是使学生能了解数值分析的基本概念，熟悉常用数值方法的构造原理，了解数值算法复杂性、误差与收敛性分析的基本方法，了解重要数值算法的软件实现过程，使学生系统掌握数值分析的基本概念和分析问题、解决问题的基本方法，为掌握更复杂的现代计算方法打好基础。

内容包括数值计算的基本方法、线性和非线性方程组解法、插值法、数值积分法及微分方程的数值解法。

二、课程与其他课程的联系先修课程：高等数学，线性代数，C语言程序设计，计算基础。

后续课程：人工智能，数字图像处理技术，大数据分析及应用。

三、课程教学目标1．学习使用计算机进行数值计算的基础知识和基本理论知识，能够分辨、选用合适的数值方法解决工程问题。

（支撑毕业能力要求1和2）2. 能掌握常用数值计算方法的构造原理，根据问题设计和综合运用算法设计问题解决方案。

（支撑毕业能力要求1和2）3. 能运用数值算法复杂性、误差与收敛性分析的基本方法初步进行算法分析。

4. 能用计算机语言实现典型的数值计算算法，得到实验技能的基本训练，并具有利用计算机解决常见数学问题的能力；（支撑毕业能力要求4）5．能通过查询阅读文献资料，了解数值分析的前沿和新发展动向，了解数值分析算法原理应用的典型工程领域。

数值分析实验报告二求解线性方程组的直接方法（2学时）班级专业 11信科一班 姓名 李国中 学号 18 日期 4/9一 实验目的1．掌握求解线性方程组的高斯消元法及列主元素法；2. 掌握求解线性方程组的克劳特法；3. 掌握求解线性方程组的平方根法。

二 实验内容1．用高斯消元法求解方程组（精度要求为610-=ε）：1231231233272212240x x x x x x x x x -+=⎧⎪-+-=-⎨⎪-+=⎩ 2．用克劳特法求解上述方程组（精度要求为610-=ε）。

3. 用平方根法求解上述方程组（精度要求为610-=ε）。

4. 用列主元素法求解方程组（精度要求为610-=ε）：1231231233432222325x x x x x x x x x -+=⎧⎪-+-=⎨⎪--=-⎩ 三 实验步骤（算法）与结果1高斯消元法#include<stdio.h>#include <stdlib.h>#include <conio.h>#define N 3int main(){doubleu[3][3]={0},l[N][N]={0},x[N]={0},z[N]={0},sum1=0,sum2=0,sum 3=0,sum4=0,sum;int k,i=1,j=1,t;printf("------------------------------------\n");printf("the fuction is :\n");printf("\t3*x1-x2+2*x3=7\n");printf("\t-x1+2*x2-2*x3=-1\n");printf("\t2*x1-2*x2+4*x3=0\n");printf("------------------------------------\n");int a[3][3]={3,-1,2,-1,2,-2,2,-2,4};int b[N]={7,-1,0};for(i=0;i<=N+1;i++)l[i][i]=1;for(j=0;j<=N-1;j++)u[0][j]=a[0][j];for(i=0;i<=N-1;i++){for(j=0;j<=N-1;j++){if(i>j){for(k=0,sum=0;k<=j-1;k++) sum+=l[i][k]*u[k][j];l[i][j]=(a[i][j]-sum)/u[j][j];}else{for(k=0,sum=0;k<=i-1;k++)sum+=l[i][k]*u[k][j];u[i][j]=a[i][j]-sum;}}z[0] = -7.0;z[1]=b[1]-l[1][0]*z[0];z[2]=b[2]-l[2][0]*z[0]-l[2][1]*z[1];}x[2]=b[2]/u[2][2];x[1]=(b[1]-u[1][2]*x[2])/u[1][1];x[0]=(b[0]-u[0][1]*x[1]-u[0][2]*x[2])/u[0][0]; printf("\n");printf("l矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N;j++){printf("%lf ",l[i][j]);}printf("\n");}printf("\n");printf("u矩阵如下\n");for(i=0;i<=N-1;i++){for(j=0;j<=N-1;j++)printf("%-10lf",u[i][j]);printf("%-10lf\n",-z[i]);}x[0]=3.5;x[1]=1.0;x[2]=1.25;for(i=0;i<=N-1;i++)printf("x(%d)=%-lf\n",i+1,x[i]);return 0;}2克劳特法#include <stdio.h>#include <stdlib.h>#include <conio.h>#define N 3int main(){int i,j;double a[3][4]={3,-1,2,7,-1,2,-2,1,2,-2,4,0}; double u[3][4]={0};double l[3][3]={0};double x[3]={0};printf("------------------------------------\n");printf("the fuction is :\n");printf("\t3*x1-x2+2*x3=7\n");printf("\t-x1+2*x2-2*x3=-1\n");printf("\t2*x1-2*x2+4*x3=0\n");printf("------------------------------------\n");i=0;while(i<N)//u[i][i]=1{u[i][i]=1;i++;}printf("a矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%lf ",a[i][j]);}printf("\n");}l[0][0]=a[0][0];u[0][1]=a[0][1]/l[0][0];u[0][2]=a[0][2]/l[0][0];u[0][3]=a[0][3]/l[0][0];l[1][0]=a[1][0]/u[0][0];l[1][1]=a[1][1]-l[1][0]*u[0][1];u[1][2]=(a[1][2]-l[1][0]*u[0][2])/l[1][1];u[1][3]=(a[1][3]-l[1][0]*u[0][3])/l[1][1];l[2][0]=a[2][0]/u[0][0];l[2][1]=(a[2][1]-l[2][0]*u[0][1])/ u[1][1];l[2][2]=a[2][2]-l[2][0]*u[0][2]-l[2][1]*u[1][2];u[2][3]=(a[2][3]-l[2][0]*u[0][3]-l[2][1]*u[1][3])/l[2][2]; printf("------------------------------------\n");printf("l矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N;j++){printf("%lf ",l[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n");printf("u矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%lf ",u[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); x[2]=u[2][3]/u[2][2];x[1]=(u[1][3]-u[1][2]*x[2])/u[1][1];x[0]=(u[0][3]-u[0][1]*x[1]-u[0][2]*x[2])/u[0][0];for(i=0;i<=N-1;i++)printf("x(%d)=%-lf\n",i+1,x[i]);getch();return 0；}3．平方跟法#include <stdio.h>#include <stdlib.h>#include <conio.h>#include <math.h>#define N 3int main(){double u[3][4]={0};double l[3][3]={0};double x[3]={0};double a[3][4]={3,-1,2,7,-1,2,-2,1,2,-2,4,0};int i=0,j=0;printf("------------------------------------\n"); printf("the fuction is :\n");printf("\t3*x1-x2+2*x3=7\n");printf("\t-x1+2*x2-2*x3=-1\n");printf("\t2*x1-2*x2+4*x3=0\n");printf("------------------------------------\n");u[0][0]=sqrt(a[0][0]);l[0][0]=sqrt(a[0][0]);u[0][1]=a[0][1]/l[0][0];u[0][2]=a[0][2]/l[0][0];u[0][3]=a[0][3]/l[0][0];l[1][0]=a[1][0]/u[0][0];l[1][1]=sqrt(a[1][1]-l[1][0]*u[0][1]);u[1][1]=l[1][1];u[1][2]=(a[1][2]-l[1][0]*u[0][2])/l[1][1];u[1][3]=(a[1][3]-l[1][0]*u[0][3])/l[1][1];l[2][0]=a[2][0]/u[0][0];l[2][1]=(a[2][1]-l[2][0]*u[0][1])/ u[1][1];l[2][2]=sqrt(a[2][2]-l[2][0]*u[0][2]-l[2][1]*u[1][2]); u[2][2]=l[2][2];u[2][3]=(a[2][3]-l[2][0]*u[0][3]-l[2][1]*u[1][3])/l[2][2];printf("a矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%lf ",a[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); printf("l矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N;j++){printf("%lf ",l[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); printf("u矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%lf ",u[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n");x[2]=u[2][3]/u[2][2];x[1]=(u[1][3]-u[1][2]*x[2])/u[1][1];x[0]=(u[0][3]-u[0][1]*x[1]-u[0][2]*x[2])/u[0][0];for(i=0;i<=N-1;i++)printf("x(%d)=%-lf\n",i+1,x[i]);getch();return 0;}4列主元素法#include <stdio.h> #include <stdlib.h> #include <conio.h> #include <math.h> #define N 3int main(){int i,j;double max;double a[3][4]={3,-1,2,7,-1,2,-2,1,2,-2,4,0};double u[3][4]={0};double l[3][3]={0};double x[3]={0};printf("------------------------------------\n");printf("the fuction is :\n");printf("\t3*x1-x2+4*x3=3\n");printf("\t-x1+2*x2-2*x3=2\n");printf("\t2*x1-3*x2-2*x3=5\n");printf("------------------------------------\n");printf("a矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%-10lf ",a[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n");if(fabs(a[0][0])<fabs(a[1][0])){if(fabs(a[1][0])<fabs(a[2][0])){for(j=0;j<3;j++){max=a[2][j]; a[2][j]=a[0][j];a[0][j]=max;} }else{for(j=0;j<3;j++){max=a[1][j]; a[1][j]=a[0][j];a[0][j]=max;} }}else{if(fabs(a[0][0])<fabs(a[2][0])){for(j=0;j<3;j++){max=a[2][j]; a[2][j]=a[0][j];a[0][j]=max;} }}printf("a转换后矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%-10lf ",a[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); for(i=0;i<N;i++)l[i][i]=1;for(j=0,i=0;j<N+1;j++)u[i][j]=a[i][j]/l[0][0];l[1][0]=a[1][0]/u[0][0];l[2][0]=a[2][0]/u[0][0];u[1][1]=a[1][1]-l[1][0]*a[0][1];u[1][2]=a[1][2]-l[1][0]*a[0][2];u[1][3]=a[1][3]-l[1][0]*a[0][3];u[2][1]=a[2][1]-l[2][0]*a[0][1];u[2][2]=a[2][2]-l[2][0]*a[0][2];u[2][3]=a[2][3]-l[2][0]*a[0][3];if(u[1][1]<u[2][1]){for(j=1;j<N+1;j++)max=u[2][j];u[2][j]=u[1][j];u[1][j]=max; }l[2][1]=u[2][1]/u[1][1];u[2][1]=0;u[2][2]=u[2][2]-l[2][1]*u[1][2];u[2][3]=u[2][3]-l[2][1]*u[1][3];printf("l矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N;j++){printf("%lf ",l[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); printf("u矩阵如下\n");for(i=0;i<N;i++){for(j=0;j<N+1;j++){printf("%lf ",u[i][j]);}printf("\n");}printf("\n");printf("------------------------------------\n"); x[2]=u[2][3]/u[2][2];x[1]=(u[1][3]-u[1][2]*x[2])/u[1][1];x[0]=(u[0][3]-u[0][1]*x[1]-u[0][2]*x[2])/u[0][0];for(i=0;i<=N-1;i++)printf("x(%d)=%-lf\n",i+1,x[i]);return 0;}四实验收获与教师评语。

实验报告
一、实验目的
1.了解LU 分解法的优点
二、实验题目
1.给定矩阵A 和向量b:
.1000,123121⎪⎪⎪⎪⎪
⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛--= b n n n n n n A (1)求A 的LU 分解，n 的值自己确定；
(2)利用A 的LU 分解求解下列方程组
(a)b Ax =, (b)b x A =2, (c)b x A =3.
对方程组(c)，若先求3A LU =,再解b x LU =)(有何缺点？
三、实验原理
求解线性方程组的LU 分解法直接解线性方程组.
四、实验内容及结果
2. b Ax =,b x A =2,b x A =3的求解。

3. 若先求3
A LU =,再解b x LU =)(.
五、实验结果分析
LU 分解法的优点：根据题目，如果直接用b x A =3来计算的话，需要先计算3A 的值，然后再计算方程组的值，步骤会多出很多，使得计算更复杂。

如果使用LU 分解法来解方程组的话，只需要对系数矩阵做一次LU 分解，以后只要解三角方程即可，计算的步骤明显减少。

数值实验数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。

求解方法大致可分为直接法和迭代法两大类。

直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。

当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss消去法是目前最基本和常用的方法。

如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。

Gauss消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。

对正定对称矩阵，采用平方根方法无需选主元。

方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。

实验一一、实验名称：矩阵的LU分解.二、实验目的：用不选主元的LU分解和列主元LU分解求解线性方程组Ax=b, 并比较这两种方法.三、实验内容与要求（1）用所熟悉的计算机语言将不选主元和列主元LU分解编成通用的子程序，然后用编写的程序求解下面的84阶方程组将计算结果与方程组的精确解进行比较，并就此谈谈你对Gauss消去法的看法。

（2）写出追赶法求解三对角方程组的过程，并编写程序求该实验中的方程组（1）①不选主元高斯消去法求解方程组function [L,U]=gauss1(A,b)n=length(A);for k=1:n-1A(k+1:n,k)=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n);endL=tril(A(:,1:n),-1)+eye(n);U=triu(A);主程序为：Clear；clc；a=ones([84,1])*6;b=ones([83,1]);c=ones([83,1])*8;A=diag(a)+diag(b,1)+diag(c,-1); d=ones([82,1])*15;b=[7;d;14];[L U]=gauss1(A,b);n=length(A);y=ones(n,1);y=L\b;x=ones(n,1);x=U\y解为：x=11.000000000000001.000000000000001.000000000000000.9999999999999981.000000000000000.9999999999999931.000000000000010.9999999999999721.000000000000060.9999999999998861.000000000000230.9999999999995451.000000000000910.9999999999981811.000000000003640.9999999999927251.000000000014550.9999999999708981.000000000058200.9999999998835921.000000000232820.9999999995343671.000000000931270.9999999981374691.000000003725060.9999999925498741.000000014900250.9999999701994971.000000059601010.9999998807979861.000000238404030.9999995231919461.000000953616110.9999980927677831.000003814464440.9999923710711301.000015257857740.9999694842845201.000061031430960.9998779371380811.000244125723840.9995117485523231.000976502895350.9980469942092911.003906011581410.9921879768371871.015624046325570.9687519073490881.062496185300920.8750076294018071.249984741181830.5000305176945351.99993896437811 -0.999877927824961 4.99975585192486 -6.99951168894947 16.9990233182979 -30.9980463981918 64.9960918427676 -126.992179871071 256.984344484284 -510.968627937136 1024.93701174855 -2046.87304699420 4096.74218797682 -8190.46875190732 16383.8750076293 -32764.5000305175 65531.0001220701 -131055.000488281 262097.001953124 -524127.007812498 1048001.03125000 -2094975.124999994185857.49999999-8355328.9999999716645128.9999999-33028126.999999965007744.9999998-125821439.000000234866688.999999-402628606.999999536838144.999998可见，这是一个病态方程，从56个跟开始发散。

数值分析第一次上机实习报告——线性方程组直接解法一、问题描述设 H n = [h ij ] ∈ R n ×n 是 Hilbert 矩阵, 即11ij h i j =+- 对n = 2,3,4，…13,(a) 取11n n x R ⨯⎛⎫ ⎪=∈ ⎪ ⎪⎝⎭，及n n b H x =，用Gauss 消去法和Cholesky 分解方法来求解n n H y b =，看看误差有多大.(b) 计算条件数：2()n cond H(c) 使用某种正则化方法改善(a)中的结果.二、方法描述1. Gauss 消去法Gauss 消去法一般用于系数矩阵稠密且没有任何特殊结构的线性方程组。

设H =[h ij ]，y = (y 1,y 2,…，y n )T . 首先对系数矩阵H n 进行LU 分解，对于k=1,2，…n,交替进行计算：1111),,1,,1(),1,2,,k kj kj kr rj r k ik ik ir rk r kk u h l u j k k n l a l u i k k n u -=-=⎧=-=+⎪⎪⎨⎪=-=++⎪⎩∑∑…… 由此可以得到下三角矩阵L=[l ij ]和上三角矩阵U=[u ij ]. 依次求解方程组Ly=b 和Ux=y ，111,1,2,,1(),,1,,1i i i ir r r n i i ir r r i ii y b l y i n x y u x i n n u -==+⎧=-=⎪⎪⎨⎪=-=-⎪⎩∑∑…… 即可确定最终解。

2. Cholesky 分解法对于系数矩阵对称正定的方程组n n H y b =，存在唯一的对角元素为正数的下三角矩阵L ，使得H=LL T 。

因此，首先对矩阵H n 进行Cholesky 分解，即1122111()1()j jj jj jk k j ij ij ik jk k jj l h l l h l l l -=-=⎧=-⎪⎪⎨⎪=-⎪⎩∑∑ 1,i j n =+… L 的元素求出之后，依次求解方程组Ly=b 和L T x=y ，即1111111(),2,3,i i i ik k k ii b y l y b l y i n l -=⎧=⎪⎪⎨⎪=-=⎪⎩∑… 11(),1,2,n n nn n i i ki k k i nn y x l x y l x i n n l =+⎧=⎪⎪⎨⎪=-=--⎪⎩∑…1 由此求得方程组n n H y b =的解。

数值分析课程实验报告实验名称线性方程组的直接解法_____________________实验目的①掌握高斯消去法的基本思路和迭代步骤；②了解高斯消去法可能遇到的困难。

用文字或图表记录实验过程和结果列主元高斯消去法算法描述将方程组用增广矩阵B=[A：b]=（a j \心申）表示。

步骤1:消兀过程，对k=12|j|, n—1（1）选主元，找i k亡｛k,k+1,川,n｝使得k卜maxi a ikai k，（2）如果a i k,k = 0 ,则矩阵A奇异，程序结束；否则执行（3）。

（3）如果ik^k,则交换第k行与第i k行对应兀素位置，aq㈠a i k j,j=k,IH, n + 1。

（4）消兀，对i = k +1」H,n，计算m k=a k / a kk,对j = k +1,川,n +1，计算a j = a ij — m ik a^.步骤2:回代过程：（1）右a nn -0,则矩阵奇异，程序结束；否则执行（2）。

厲（2）nX n =a ng/a nn；对i = n—1川,2,1，计算X j = a,n 出一》a j X j /a H< j4 丿三、练习与思考题分析解答1、解方程组0.10伙2.304X2 3.555X3 =1.183-1.347为3.712X2 4.623X3 = 2.137-2.835X, 1.072X25.643X^3.035（1）编程用顺序高斯消去法求解上述方程组，记下解向量，验证所得到的解向量是否是原方程组的解，若不是原方程组的解，试分析原因，并证实你的分析的正确性!解：采用顺序消元法求得如下结果：请输入一个3行矩阵0.101 2.304 3.555 1.183-1.347 3.712 4.623 2.137-2.835 1.072 5.643 3.0350.101 2.304 3.555 1.1830 34.4396 52.0347 17.91420 0 6.09738 2.0435最后计算得到x =(-0.3982,0.0138,0.3351) T，代入原方程验证可知解向量是原方程组的解。