数值分析实验2_求解线性方程组直接法
- 格式:doc
- 大小:58.50 KB
- 文档页数:6
《数值分析》课程教学大纲课程编号: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,代入原方程验证可知解向量是原方程组的解。
广东金融学院实验报告课程名称:数值分析实验目的及要求实验目的:题一:通过数值实验,从中体会解线性方程组选主元的必要性和LU分解法的优点,以及方程组系数矩阵和右端向最的微小变化对解向最的影响。
比较各种直接接法在解线性方程组中的效果;题二:认识齐种迭代法收敛的含义、影响齐迭代法收敛速度的因素。
实验要求:题一:(1)在MATLAB中编写程序用列主元高斯消去法和LU分解求解上述方程组,输出曲b中矩阵A 及向量b和A二LU分解中的L及U, detA及解向量X.(2)将方程组中的2. 099999改为2. 1, 5. 900001改为5. 9,用列主元高斯消去法求解变换后的方程组,输出解向最x及detA,并与(1)中的结果比较。
(3)用MATLAB的内部函数inv求出系数矩阵的逆矩阵,再输入命令x=inv(A)*b,即可求出方程组的解。
请与列主元高斯消公法和LU分解法求出的解进行比较,体会选主元的方法具有良好的数值稳定性。
用MATLAB的内部曲数det求出系数行列式的值,并与(1)、(2)中输出的系数行列式的值进行比较。
(4)比较以上各种直接解法在解线性方程组中的效果。
题二:(1)选取不同的初始向M:X(0)及右端向最b,给泄迭代误差要求,用Jacobi迭代法和Gauss-Seidel迭代法求解,观察得到的序列是否收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论。
列岀算法清单。
(2)用SOR迭代法求上述方程组的解,松弛系数血取1<69<2的不同的三个值,在< 10"5时停止迭代,记录迭代次数,分析计算结呆与松弛系数血的关系并得出你的结论。
(3)用MATLAB的内部函数inv求出系数矩阵的逆矩阵.再输入命令^inv(A)*b>即可求出上述各个方程组的解.并与上述三种方法求出的解进行比较。
请将比较结果列入卜表。
方程组的解X1 Xr■迭代次数误差楮确解Jacibi解法Gause・seidel 解法SOR 解法00= 60= 60=实验环境及相关情况(包含使用软件、实验设备、主要仪器及材料等)1. Win72. Mat lab 7.0实验内容及步骤(包含简要的实验步骤流程) 实验内容:题一:解卜列线性方程组'10 -7‘X 】、(8、-3 2.099999 62Xr5.9000015-1 5 -15、12> 0< 1 >题二研究解线性方程组 做=b 迭代法的收敛性、收敛速度以及SOR 方法中/佳松弛因子的选取问题, 用迭代法求解做二b,其中・4 -1r■7 A=4 -81 ,b =-21-2 ■1515实验结果(包括程序或图表、结论陈述.数据记录及分析等,可附页)题一:直接解法解线性方程组(1)列主兀高斯消去法与LU 分解求解列主元高斯消去法:编写matalab 程序(见附录gaosi.m ),输出矩阵10.000 -7.000 0.000= 0.000 2.5000-5.000一 0.000 0.0006.0000020.000 0.000 0.000向量8 b =1 8.300 L5.0800J解向量:X = (0 ・-1 , 1 r I )7 其中系数行列式的值det (A )=762.00009LU 分解求解:编J matalab 程序(见附录zhjLU. m 和LU ・m ),执行输出:-1.5 2.300 5.080-3.0001.000000.00000.5000 -25000001.0000 0.2000 -24000000.9600 10.0000 -7.0000 0.0000 1.0000n = 0.0000-0.0000010.0000 2.3000 —0.0000 0.000015000000 57500000.0000 0.0000 0.0000 5.0800在matlab 命令窗II 输入L*U ,可以得到A 二L*U ,即分解结果正确。
数值分析实验报告实验一、解线性方程组的直接方法——梯形电阻电路问题利用追赶法求解三对角方程组的方法,解决梯形电阻电路问题:电路中的各个电流{1i ,2i ,…,8i }须满足下列线性方程组:R V i i =- 22 210 252321=-+-i i i 0 252 432=-+-i i i 0 252 543=-+-i i i 0 252 654=-+-i i i 0 252 765=-+-i i i 0 252 876=-+-i i i 052 87=+-i i设V 220=V ,Ω=27R ,运用追赶法,求各段电路的电流量。
问题分析:上述方程组可用矩阵表示为:⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡--------------00000001481.8522520000002520000002520000002520000002520000002520000002287654321i i i i i i i i问题转化为求解A x b =,8阶方阵A 满足顺序主子式(1,2...7)0i A i =≠,因此矩阵A存在唯一的Doolittle 分解,可以采用解三对角矩阵的追赶法!追赶法a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0]; d=[220/27 0 0 0 0 0 0 0];Matlab 程序function x= zhuiganfa( a,b,c,d )%追赶法实现要求:|b1|>|C1|>0,|bi|>=|ai|+|ci| n=length(b); u=ones(1,n); L=ones(1,n); y=ones(1,n); u(1)=b(1); y(1)=d(1); for i=2:nL(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*L(i); y(i)=d(i)-y(i-1)*L(i); endx(n)=y(n)/u(n); for k=n-1:-1:1x(k)=(y(k)-c(k)*x(k+1))/u(k); end endMATLAB 命令窗口输入:a=[0 -2 -2 -2 -2 -2 -2 -2]; b=[2 5 5 5 5 5 5 5];c=[-2 -2 -2 -2 -2 -2 -2 0] d=[220/27 0 0 0 0 0 0 0];x= zhuiganfa(a,b,c,d )运行结果为:x =8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477存在问题根据电路分析中的所讲到的回路电流法,可以列出8个以回路电流为独立变量的方程,课本上给出的第八个回路电流方程存在问题,正确的应该是78240i i -+=;或者可以根据电路并联分流的知识,同样可以确定78240i i -+=。
《数值分析》实验报告班级:姓名:学号:指导老师:实验基本要求一、上机前的准备工作1、复习和掌握与本次实验有关的教学内容。
2、根据本次实验要求,在纸上编写算法及上机的程序,并经过人工模拟运行检验,减少不必要的错误,提高上机效率。
切忌不编程序、不作人工检查就进行程序输入,这只能使上机调试的难度增加,甚至可能带来学习自信心的下降,影响后续课程的学习。
二、上机实验步骤1、启动开发环境;2、建立源程序文件,输入源程序;3、编译产生目标程序,连接生成可执行程序,运行程序,输出结果;4、对数值计算结果进行误差分析,讨论数值算法的收敛性与稳定性;5、整理实验报告。
三、实验报告实验报告是记录实验工作全过程的技术文档,实验报告的撰写是科学技术工作的一个组成部分。
《数值分析》实验报告包括下列要求:1、实验原理;2、实验内容和要求;3、数值算法描述,包括数据输入、数据处理和数据输出;4、算法的实现(1)给出具体的计算实例,(2)经调试正确的源程序清单,(3)对具体的数值例子给出数值结果;5、计算结果的误差分析,算法的收敛性与稳定性的讨论;6、实验心得。
实验一、误差分析一、实验目的1、通过上机编程,复习巩固以前所学程序设计语言及上机操作指令;2、通过上机计算,了解误差、绝对误差、误差界、相对误差界的有关概念;3、 通过上机计算,了解舍入误差所引起的数值不稳定性。
二、实验原理误差问题是数值分析的基础,又是数值分析中一个困难的课题。
在实际计算中,如果选用了不同的算法,由于舍入误差的影响,将会得到截然不同的结果。
因此,选取算法时注重分析舍入误差的影响,在实际计算中是十分重要的。
同时,由于在数值求解过程中用有限的过程代替无限的过程会产生截断误差,因此算法的好坏会影响到数值结果的精度。
三、实验任务对20,,2,1,0 =n ,计算定积分⎰+=105dx x x y nn . 算法1:利用递推公式151--=n n y ny , 20,,2,1 =n , 取 ⎰≈-=+=100182322.05ln 6ln 51dx x y . 算法2:利用递推公式n n y n y 51511-=- 1,,19,20 =n . 注意到⎰⎰⎰=≤+≤=1010202010201051515611261dx x dx x x dx x , 取 008730.0)12611051(20120≈+≈y . 思考:从计算结果看,哪个算法是不稳定的,哪个算法是稳定的。
数值分析实验报告二求解线性方程组的直接方法(2学时)一 实验目的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. 高斯消元法求解根据算法思想用C 语言编程(源程序见附录2.1) 编译结果如下图:2. 克劳特法求解根据算法思想用C 语言编程(源程序见附录2.2)编译结果如下图:3平方根法求解根据算法思想用C语言编程(源程序见附录2.3)编译结果如下图:4列主元素法求解根据算法思想用C语言编程(源程序见附录2.4)编译结果如下图:四实验收获与教师评语1.实验收获:对于这次实验,我可以锻炼到上机实验的能力,并且第一次感受到数学知识在现实生活中的应用,也是第一次运用计算机解决数学问题。
另外,正是因为这次上机实验,让我重温了有些遗忘的编程知识。
2.教师评语:附录:2.1至2.4源程序代码2.1 #include <stdio.h>#define N 8void main(){floatsum,a[N][N]={0},u[N][N]={0},l[N][ N]={0},z[N]={0},x[N]={0};int n,i,j,k;printf("input the number of roots:");/*方程的个数小于8*/scanf("%d",&n);printf("input xishu matrix:\n"); for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)scanf("%f",&a[i][j]);for(i=1;i<=n;i++)l[i][i]=1;for(i=1;i<=n;i++)for(j=1;j<=n;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; }}for(i=1,j=n+1;i<=n;i++){for(k=0,sum=0;k<=i;k++)sum+=l[i][k]*z[k];z[i]=a[i][j]-sum;}for(i=n;i>=1;i--){sum=0;for(k=i+1;k<=n;k++)sum+=u[i][k]*x[k]; x[i]=(z[i]-sum)/u[i][i];}printf("changing fangcheng xi shu is:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("%-10g",u[i][j]); printf("%-10g\n",z[i]);}for(i=1;i<=n;i++)printf("x(%d)=%-g\n",i,x[i]); getch();}2.2#include <conio.h>#include <stdio.h>#define N 8void main(){floatsum,a[N][N]={0},u[N][N]={0},l[N][ N]={0},z[N]={0},x[N]={0};int n,i,j,k;printf("input the number of roots:");/*方程的个数小于8*/scanf("%d",&n);printf("input xishu matrix:\n"); for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)scanf("%f",&a[i][j]);for(i=1;i<=n;i++)u[i][i]=1;for(i=1;i<=n;i++)for(j=1;j<=n;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; }else{for(k=0,sum=0;k<=i-1;k++)sum+=l[i][k]*u[k][j];u[i][j]=(a[i][j]-sum)/l[i][i];}}for(i=1,j=n+1;i<=n;i++){for(k=0,sum=0;k<=i;k++)sum+=l[i][k]*z[k];z[i]=(a[i][j]-sum)/l[i][i];}for(i=n;i>=1;i--){sum=0;for(k=i+1;k<=n;k++)sum+=u[i][k]*x[k];x[i]=z[i]-sum;}printf("changing fangcheng xi shu is:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("%-10g",u[i][j]); printf("%-10g\n",z[i]);}for(i=1;i<=n;i++)printf("x(%d)=%-g\n",i,x[i]); getch();}2.3#include <conio.h>#include <stdio.h>#include <math.h>#define N 8void main(){floatsum,a[N][N]={0},u[N][N]={0},l[N][ N]={0},z[N]={0},x[N]={0};int n,i,j,k;printf("input the number of roots:");/*方程的个数小于8*/ scanf("%d",&n);printf("input xishu matrix:\n"); for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)scanf("%f",&a[i][j]);for(i=1;i<=n;i++)for(j=1;j<=n;j++){if(i>j){for(k=0,sum=0;k<=j-1;k++)sum+=l[i][k]*l[j][k];l[i][j]=(a[i][j]-sum)/l[j][j];u[j][i]=l[i][j];}if(i==j){for(k=0,sum=0;k<=i-1;k++)sum+=l[i][k]*l[i][k];l[i][i]=sqrt(a[i][i]-sum);u[i][i]=l[i][i];}}for(i=1,j=n+1;i<=n;i++){for(k=0,sum=0;k<=i;k++)sum+=l[i][k]*z[k];z[i]=(a[i][j]-sum)/l[i][i];}for(i=n;i>=1;i--){sum=0;for(k=i+1;k<=n;k++)sum+=u[i][k]*x[k];x[i]=(z[i]-sum)/l[i][i];}printf("changing fangcheng xi shu is:\n");for(i=1;i<=n;i++){for(j=1;j<=n;j++)printf("%-10g",u[i][j]); printf("%-10g\n",z[i]);}for(i=1;i<=n;i++)printf("x(%d)=%-g\n",i,x[i]); getch();}2.4#include <stdio.h>#include <math.h>#define N 8void main(){ int i,j,k,n;floata[N][N],l[N][N],u[N][N],z[N],x1[N ],max,x;for(i=0;i<N;i++)for(j=0;j<N;j++)a[i][j]=u[i][j]=l[i][j]=0;z[0]=0; printf("input the number of roots:");scanf("%d",&n);printf("input xishu matrix:\n");for(i=1;i<=n;i++)for(j=1;j<=n+1;j++)scanf("%f",&a[i][j]);for(i=1;i<n;i++){ max=fabs(a[i][i]);k=i;for(j=i+1;j<=n;j++)if(fabs(a[j][i])>max){ max=fabs(a[j][i]);k=j;} for(j=i;j<=n+1;j++){ x=a[i][j];a[i][j]=a[k][j];a[k][j]=x;}for(j=i+1;j<=n;j++){ l[j][i]=a[j][i]/a[i][i];a[j][i]=0;for(k=i+1;k<=n+1;k++) a[j][k]=a[j][k]-l[j][i]*a[i][k]; }z[i]=a[i][n+1];}z[n]=a[n][n+1];for(i=1;i<=n;i++)for(j=1;j<=n;j++)u[i][j]=a[i][j];for(i=n;i>0;i--){ k=i+1;x=0;while(k<=n){ x+=u[i][k]*x1[k];k++;}x1[i]=(z[i]-x)/u[i][i];}for(i=1;i<=n;i++)for(j=1;j<=n;j++){ printf("%-10g",u[i][j]);if(j==n)printf("%-10g\n",z[i]);}for(i=1;i<=n;i++)printf("x%d=%-10g\n",i,x1[i]);getch();}。
数值分析第五章解线性方程组的直接法解线性方程组是数值分析中的一个重要问题,对于大规模的线性方程组来说,直接法是一种常用的求解方法。
本文将介绍解线性方程组的直接法,包括高斯消元法和LU分解法,并对其稳定性和计算复杂度进行讨论。
高斯消元法是一种常用的直接法,用于求解非奇异线性方程组。
其基本思想是通过初等行变换将线性方程组转化为上三角方程组,然后通过回代求解得到方程的解。
高斯消元法的步骤如下:1.将线性方程组表示为增广矩阵[A,b],其中A是系数矩阵,b是常数向量。
2.从第一行开始,选择一个非零元素作为主元,通过行变换将主元下方的元素全部消为零。
3.重复第2步,直到矩阵变为上三角矩阵。
4.通过回代求解上三角矩阵,得到方程组的解。
高斯消元法的主要优点是简单直接,容易实现,但存在一些问题。
首先,如果系数矩阵A是奇异矩阵,即行列式为零,那么高斯消元法无法得到方程组的解。
其次,如果系数矩阵A的其中一行或几行接近于线性相关,那么在消元过程中会引入大量的舍入误差,导致计算结果不准确。
这也说明了高斯消元法的稳定性较差。
为了提高稳定性,可以使用LU分解法来解线性方程组。
LU分解法将系数矩阵A分解为两个矩阵L和U的乘积,其中L是下三角矩阵,U是上三角矩阵。
这样,原始的线性方程组可以表示为LUx=b,进而可以通过两个步骤来求解方程组:1.进行LU分解,将系数矩阵A分解为L和U。
2.分别用前代和回代的方法求解方程组Ly=b和Ux=y。
LU分解法相对于高斯消元法的优点是,可以在求解多个右端向量时,避免重复计算LU分解,从而提高计算效率。
同时,LU分解法的稳定性也较高,对于多个右端向量求解时,舍入误差的累积相对较小。
然而,LU分解法也存在一些问题。
首先,LU分解法的计算复杂度较高,需要进行两次矩阵乘法和一次矩阵向量乘法,而且LU分解过程中需要对系数矩阵A进行大量的行变换,增加了计算量。
其次,当系数矩阵A的一些元素非常小或非常大时,LU分解法容易出现数值不稳定的情况,即舍入误差的累积较大,导致计算结果不准确。
实验五 解线性方程组的直接方法实验5.1 (主元的选取与算法的稳定性) 问题提出:Gauss 消去法是我们在线性代数中已经熟悉的。
但由于计算机的数值运算是在一个有限的浮点数集合上进行的,如何才能确保Gauss 消去法作为数值算法的稳定性呢?Gauss 消去法从理论算法到数值算法,其关键是主元的选择。
主元的选择从数学理论上看起来平凡,它却是数值分析中十分典型的问题。
实验内容:考虑线性方程组n n n R b R A b Ax ∈∈=⨯,,编制一个能自动选取主元,又能手动选取主元的求解线性方程组的Gauss 消去过程。
实验要求:(1)取矩阵⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=1415157,6816816816 b A ,则方程有解T x )1,,1,1(* =。
取n=10计算矩阵的条件数。
让程序自动选取主元,结果如何?(2)现选择程序中手动选取主元的功能。
每步消去过程总选取按模最小或按模尽可能小的元素作为主元,观察并记录计算结果。
若每步消去过程总选取按模最大的元素作为主元,结果又如何?分析实验的结果。
(3)取矩阵阶数n=20或者更大,重复上述实验过程,观察记录并分析不同的问题及消去过程中选择不同的主元时计算结果的差异,说明主元素的选取在消去过程中的作用。
(4)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
思考题一:(Vadermonde 矩阵)设⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=∑∑∑∑====n i i n n i i ni i n i i n n n n n n nx x x x b x x x x x x x x x x x x A 002010022222121102001111 ,, 其中,n k k x k ,,1,0,1.01 =+=,(1)对n=2,5,8,计算A 的条件数;随n 增大,矩阵性态如何变化?(2)对n=5,解方程组Ax=b ;设A 的最后一个元素有扰动10-4,再求解Ax=b(3)计算(2)扰动相对误差与解的相对偏差,分析它们与条件数的关系。
《数值分析》课程实验报告范文《数值分析》课程实验报告姓名:学号:学院:机电学院日期:2022年某月某日目录实验一函数插值方法1实验二函数逼近与曲线拟合5实验三数值积分与数值微分7实验四线方程组的直接解法9实验五解线性方程组的迭代法15实验六非线性方程求根19实验七矩阵特征值问题计算21实验八常微分方程初值问题数值解法24实验一函数插值方法一、问题提出对于给定的一元函数的n+1个节点值。
试用Lagrange公式求其插值多项式或分段二次Lagrange插值多项式。
实验二函数逼近与曲线拟合一、问题提出从随机的数据中找出其规律性,给出其近似表达式的问题,在生产实践和科学实验中大量存在,通常利用数据的最小二乘法求得拟合曲线。
在某冶炼过程中,根据统计数据的含碳量与时间关系,试求含碳量与时间t的拟合曲线。
t(分)051015202530354045505501.272.162.863.443.874.154.374.51 4.584.024.64二、要求1、用最小二乘法进行曲线拟合;2、近似解析表达式为;3、打印出拟合函数,并打印出与的误差,;4、另外选取一个近似表达式,尝试拟合效果的比较;5、某绘制出曲线拟合图。
三、目的和意义1、掌握曲线拟合的最小二乘法;2、最小二乘法亦可用于解超定线代数方程组;3、探索拟合函数的选择与拟合精度间的关系四、实验步骤:第一步先写出线性最小二乘法的M文件functionc=lpoly(某,y,m)n=length(某);b=zero(1:m+1);f=zero(n,m+1); fork=1:m+1f(:,k)=某.^(k-1);enda=f'某f;b=f'某y';c=a\b;c=flipud(c);第二步在命令窗口输入:>>lpoly([0,5,10,15,20,25,30,35,40,45,50,55],[0,1.27,2.16,2.86,3.44,3.87,4.15,4.37,4.51,4.58,4.02,4.64],2)回车得到:an=-0.00240.20370.2305即所求的拟合曲线为y=-0.0024某2+0.2037某+0.2305在编辑窗口输入如下命令:>>某=[0,5,10,15,20,25,30,35,40,45,50,55];>>y=-0.0024某某.^2+0.2037某某+0.2305;>>plot(某,y)命令执行得到如下图五、实验结论分析复杂实验数据时,常采用分段曲线拟合方法。
一 实验目的
1.掌握求解线性方程组的高斯消元法及列主元素法;
2. 掌握求解线性方程组的克劳特法;
3. 掌握求解线性方程组的平方根法。
二 实验内容
1.用高斯消元法求解方程组(精度要求为610-=ε):
1231231
233272212240x x x x x x x x x -+=⎧⎪-+-=-⎨⎪-+=⎩ 2.用克劳特法求解上述方程组(精度要求为610-=ε)。
3. 用平方根法求解上述方程组(精度要求为610-=ε)。
4. 用列主元素法求解方程组(精度要求为610-=ε):
1231231
233432222325x x x x x x x x x -+=⎧⎪-+-=⎨⎪--=-⎩ 三 实验步骤(算法)与结果
1.
程序代码(Python3.6):
import numpy as np
def Gauss(A,b):
n=len(b)
for i in range(n-1):
if A[i,i]!=0:
for j in range(i+1,n):
m=-A[j,i]/A[i,i]
A[j,i:n]=A[j,i:n]+m*A[i,i:n]
b[j]=b[j]+m*b[i]
for k in range(n-1,-1,-1):
b[k]=(b[k]-sum(A[k,(k+1):n]*b[(k+1):n]))/A[k,k]
print(b)
运行函数:
>>> A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=np.float) >>> b=np.array([7,-1,0],dtype=np.float)
>>> x=Gauss(A,b)
输出:
结果:解得原方程的解为x1=3.5,x2=-1,x3=-2.25
2
程序代码(Python3.6):
import numpy as np
A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float)
L=np.array([[1,0,0],[0,1,0],[0,0,1]],dtype=float)
U=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float)
b=np.array([7,-1,0],dtype=float)
y=np.array([0,0,0],dtype=float)
x=np.array([0,0,0],dtype=float)
def LU(A):
n=len(A[0])
i=0
while i<n:
j=i
while j<n:
U[i,j]=A[i,j]-sum(L[i,0:i]*U[0:i,j])
j+=1
k=i+1
while k<n:
L[k,i]=(A[k,i]-sum(L[k,0:i]*U[0:i,i]))/U[i,i]
k+=1
i+=1
print('L=',L)
print('U=',U)
def solvey(L,b):
n=len(y)
y[0]=b[0]
for i in range(1,n):
y[i]=b[i]-sum(L[i,0:i]*y[0:i])
print('y=',y)
def solvex(U,y):
n=len(x)
x[n-1]=y[n-1]/U[n-1,n-1]
for i in range(n-2,-1,-1):
x[i]=(y[i]-sum(U[i,(i+1):n]*x[(i+1):n]))/U[i,i]
print('x=',x)
运行函数:
>>> LU(A)
>>> solvey(L,b)
>>> solvex(U,y)
输出:
结果:同样解得原方程的解为x1=3.5,x2=-1,x3=-2.25
3程序代码(Python3.6):
import numpy as np
A=np.array([[3,-1,2],[-1,2,-2],[2,-2,4]],dtype=float)
L=np.array([[0,0,0],[0,0,0],[0,0,0]],dtype=float)
b=np.array([7,-1,0],dtype=float)
y=np.array([0,0,0],dtype=float)
x=np.array([0,0,0],dtype=float)
def Cholesky(A):
n=len(A[0])
for k in range(n):
L[k,k]=pow(A[k,k]-sum(L[k,0:k]*L[k,0:k]),0.5)
for i in range(k+1,n):
L[i,k]=(A[i,k]-sum(L[i,0:i]*L[k,0:i]))/L[k,k]
print('L=',L)
def solvey(L,b):
n=len(y)
for i in range(n):
y[i]=(b[i]-sum(L[i,0:i]*y[0:i]))/L[i,i]
print('y=',y)
def solvex(L,y):
n=len(x)
for i in range(n-1,-1,-1):
x[i]=(y[i]-sum(L[(i+1):n,i]*x[(i+1):n]))/L[i,i]
print('x=',x)
运行函数:
>>> Cholesky(A)
>>> solvey(L,b)
>>> solvex(L,y)
输出:
结果:同样解得原方程的解为x1=3.5,x2=-1,x3=-2.25
4
程序代码(Python3.6):
import numpy as np
A=np.array([[3,-1,4],[-1,2,-2],[2,-3,-2]],dtype=float)
b=np.array([3,2,-5],dtype=float)
def Main_Gauss(A,b):
n=len(b)
for k in range(n):
A_max=0
for i in range(k,n):
if abs(A[i,k])>A_max:
A_max=abs(A[i,k])
r=i
if A_max<1e-6:
print('系数矩阵奇异,无法求解方程!')
break
if r>k:
for j in range(k,n):
s=A[k,j]
A[k,j]=A[r,j]
A[r,j]=s
t=b[k]
b[k]=b[r]
b[r]=t
for j in range(k+1,n):
A[k,j]=A[k,j]/A[k,k]
b[k]=b[k]/A[k,k]
for i in range(n):
if i!=k:
for j in range(k+1,n):
A[i,j]=A[i,j]-A[i,k]*A[k,j]
b[i]=b[i]-A[i,k]*b[k]
print('x=',b)
运行函数:
>>> Main_Gauss(A,b)
输出:
结果:解得原方程的解为x1=1,x2=2,x3=0.5
四实验收获与教师评语
实验收获:掌握了求解线性方程组的高斯消元法、列主元素法、克劳
特法和平方根法等算法流程,以及能够运用Python、MA TLAB等语言实现并解出方程组的根。