数值分析实验报告_清华大学__线性代数方程组的数值解法
数值分析计算方法实验报告
end;
end;
X=x;
disp('迭代结果:');
X
format short;
输出结果:
因为不收敛,故出现上述情况。
4.超松弛迭代法:
%SOR法求解实验1
%w=1.45
%方程组系数矩阵
clc;
A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]
b=[10,5,-2,7]'
b=[10,5,-2,7]'
[m,n]=size(A);
if m~=n
error('矩阵A的行数和列数必须相同');
return;
end
if m~=size(b)
error('b的大小必须和A的行数或A的列数相同');
return;
end
if rank(A)~=rank([A,b])
error('A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解');
3.实验环境及实验文件存档名
写出实验环境及实验文件存档名
4.实验结果及分析
输出计算结果,结果分析和小结等。
解:1.高斯列主元消去法:
%用高斯列主元消去法解实验1
%高斯列主元消元法求解线性方程组Ax=b
%A为输入矩阵系数,b为方程组右端系数
%方程组的解保存在x变量中
format long;
A=[2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13]
return;
end
c=n+1;
A(:,c)=b;
for k=1:n-1
线性方程组的数值解法实验报告
实验报告——线性方程组的数值解法姓名:班级:学号:日期:一 实践目的1. 熟悉求解线性方程组的有关理论和方法。
2. 会编列主元消去法,全主元消去法,雅克比迭代法及高斯-赛德尔迭代法的程序。
3.通过实际计算,进一步了解各种方法的优缺点,选择合适的数值方法。
4. 进一步应用数学知识,扩展数学思维,提高编程能力。
二 问题定义及题目分析1.求解线性方程是实际中常遇到的问题。
而一般求解方法有两种,一个是直接求解法,一个是迭代法。
2.直接法是就是通过有限步四则运算求的方程准确解的方法。
但实际计算中必然存在舍入误差,因此这种方法只能得到近似解。
3.迭代法是先给一个解的初始近似值,然后按一定的法则求出更准确的解,即是用某种极限过程逐步逼近准确解的方法。
4.这次实践,将采用直接法:列主元高斯消去法,全主元高斯消去法;迭代法:雅克比迭代法和高斯-赛德尔迭代法。
三 详细设计 设有线性方程组11112211n n a x a x a x b +++=21122222n n a x a x a x b +++= 1122n n nn n n a x a x a x b +++=令A= 111212122212n n n n nn a a a aa a a a a ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ x=12n x x x ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦ b= 12n b b b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦一. 上三角方程组的求解很简单。
所以若能将方程组化为上三角方程组,那就很容易得到方程的解,高斯消去法则是一种简洁高效的方法,可以将方程组化为上三角方程组,但是这种方法必须满足每一步时0kk a =才能求解,还有当kk a 绝对值很小时,作分母会引起较大的舍入误差。
因此在消元过程中应尽量选择绝对值较大的系数作为主元素。
1.列主元消去法很好的解决这一问题。
将1x 的系数1(1)k a k n ≤≤中绝对值最大者作为主元素,交换第一行和此元素所在的行,然后主元消去非主元项1x 的系数,。
数值分析实验报告——Hilbert矩阵的求解
1 / 7 数值分析课程实验报告题目:病态线性方程组的求解理论分析表明,数值求解病态线性方程组很困难。
考虑求解如下的线性方程组的求解Hx = b ,期中,期中H 是Hilbert 矩阵,()ij n n H h ´=,11ij h i j =+-,i ,j = 1,2,…,n 1.估计矩阵的2条件数和阶数的关系2.对不同的n ,取(1,1,,1,1))nx =Î,分别用Gauss 消去,Jacobi 迭代,Gauss-seidel 迭代,SOR 迭代和共轭梯度法求解,比较结果。
3.结合计算结果,试讨论病态线性方程组的求解。
解答过程1.估计矩阵的2-条件数和阶数的关系矩阵的2-条件数定义为:1222()Cond A A A-=´,将Hilbert 矩阵带入有:1222()Cond H H H -=´调用自编的Hilbert_Cond 函数对其进行计算,取阶数n = 50,可得从,可得从1阶到50阶的2-条件数,以五位有效数字输出,其中前10项见表1。
表1.前十阶Hilbert 矩阵的2-条件数阶数1 2 3 4 5 2-条件数1 19.281 524.06 1.5514e+004 4.7661e+005 阶数6 7 8 9 10 2-条件数1.4951e+007 4.7537e+008 1.5258e+010 4.9315e+011 1.6025e+013 从表1可以看出,随着阶数每递增1,Hilbert 矩阵的2-条件数都至少增加一个数量级,但难以观察出明显的相依规律。
故考虑将这些数据点绘制在以n 为横轴、Cond (H )2为纵轴的对数坐标系中(编程用Hilbert_Cond 函数同时完成了这个功能),生成结果如图1。
图1.不同阶数下Hilbert 矩阵的2-条件数分布条件数分布由图可见,当维数较小时,在y-对数坐标系中Cond (H )2与n 有良好的线性关系;但n 超过10后,线性趋势开始波动,n 超过14后更是几乎一直趋于平稳。
数值分析实验报告
一、实验目的1. 理解数值分析的基本概念和常用算法;2. 掌握数值方法在求解实际问题中的应用;3. 培养编程能力,提高对数值分析软件的使用熟练度。
二、实验内容本次实验主要涉及以下内容:1. 拉格朗日插值法;2. 牛顿插值法;3. 线性方程组的求解方法;4. 方程求根的数值方法;5. 最小二乘法曲线拟合。
三、实验步骤1. 拉格朗日插值法(1)输入数据:给定一组数据点(x1, y1)、(x2, y2)、...、(xn, yn)。
(2)计算拉格朗日插值多项式L(x)。
(3)利用L(x)计算待求点x0的函数值y0。
2. 牛顿插值法(1)输入数据:给定一组数据点(x1, y1)、(x2, y2)、...、(xn, yn)。
(2)计算牛顿插值多项式N(x)。
(3)利用N(x)计算待求点x0的函数值y0。
3. 线性方程组的求解方法(1)输入数据:给定线性方程组的系数矩阵A和常数向量b。
(2)采用高斯消元法求解线性方程组Ax=b。
4. 方程求根的数值方法(1)输入数据:给定函数f(x)和初始值x0。
(2)采用二分法求解方程f(x)=0的根。
5. 最小二乘法曲线拟合(1)输入数据:给定一组数据点(x1, y1)、(x2, y2)、...、(xn, yn)。
(2)建立线性最小二乘模型y=F(x)。
(3)利用最小二乘法求解模型参数。
四、实验结果与分析1. 拉格朗日插值法与牛顿插值法的比较通过实验,我们发现牛顿插值法的精度高于拉格朗日插值法。
这是因为牛顿插值法在计算过程中考虑了前一项的导数信息,从而提高了插值多项式的平滑性。
2. 线性方程组的求解方法高斯消元法在求解线性方程组时,计算过程较为繁琐,但稳定性较好。
在实际应用中,可根据具体问题选择合适的方法。
3. 方程求根的数值方法二分法在求解方程时,收敛速度较慢,但具有较好的稳定性。
对于初始值的选择,应尽量接近真实根。
4. 最小二乘法曲线拟合最小二乘法在拟合曲线时,误差较小,适用于数据点较多的情况。
线性代数方程组的数值解法_百度文库
线性代数方程组的数值解法【实验目的】1. 学会用MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】【题目1】通过求解线性方程组A1x=b1和A2x=b2,理解条件数的意义和方程组的性态对解的影响。
其中A1是n阶范德蒙矩阵,即⎡1x0⎢1x1⎢A1=⎢⎢⎢⎣1xn-12x0x12 2xn-1n-1⎤ x0⎥ x1n-1⎥1,...,n-1 ,xk=1+0.1k,k=0,⎥ n-1⎥ xn-1⎥⎦A2是n阶希尔伯特矩阵,b1,b2分别是A1,A2的行和。
(1)编程构造A1(A2可直接用命令产生)和b1,b2;你能预先知道方程组A1x=和A2x=。
b2的解吗?令n=5,用左除命令求解(用预先知道的解可检验程序)b1(2)令n=5,7,9,…,计算A1,A2的条件数。
为观察它们是否病态,做以下试验:b1,b2不变,A1和A2的元素A1(n,n),A2(n,n)分别加扰动ε后求解;A1和A2不变,b1,b2的分量b1(n),分析A和b的微小扰动对解的影响。
b2(n)分别加扰动ε求解。
ε取10-1010,-8,10-6。
(3)经扰动得到的解记做x~,计算误差-x~x,与用条件数估计的误差相比较。
1.1构造A1,A2和b1,b2首先令n=5,构造出A1,A2和b1,b2。
首先运行以下程序,输出A1。
运行以下程序对A1,A2求行和:由于b1,b2分别是A1,A2的行和,所以可以预知x1=运行下列程序,用左除命令对b1,b2进行求解:得到以下结果: T。
x2=(1,1, ,1)1.2 计算条件数并观察是否为病态1.不加扰动,计算条件数。
运行以下程序:由此可知,A1,A2的条件数分别是3.574∗10, 4,766∗10。
2.b1,b2不变,A1(n,n),A2(n,n)分别加扰动(1)n=5时设x11,x12,x13分别为A1添加扰动10−10,10−8,10−6后的解。
数值分析实验报告总结
一、实验背景数值分析是研究数值计算方法及其理论的学科,是计算机科学、数学、物理学等领域的重要基础。
为了提高自身对数值分析理论和方法的理解,我们进行了数值分析实验,通过实验加深对理论知识的掌握,提高实际操作能力。
二、实验目的1. 理解数值分析的基本理论和方法;2. 掌握数值分析实验的基本步骤和技巧;3. 培养实验设计和数据分析能力;4. 提高编程和计算能力。
三、实验内容本次实验主要分为以下几个部分:1. 线性方程组求解实验:通过高斯消元法、LU分解法等求解线性方程组,并分析算法的稳定性和误差;2. 矩阵特征值问题计算实验:利用幂法、逆幂法等计算矩阵的特征值和特征向量,分析算法的收敛性和精度;3. 非线性方程求根实验:运用二分法、牛顿法、不动点迭代法等求解非线性方程的根,比较不同算法的优缺点;4. 函数插值实验:运用拉格朗日插值、牛顿插值等方法对给定的函数进行插值,分析插值误差;5. 常微分方程初值问题数值解法实验:运用欧拉法、改进的欧拉法、龙格-库塔法等求解常微分方程初值问题,比较不同算法的稳定性和精度。
四、实验过程1. 线性方程组求解实验:首先,编写程序实现高斯消元法、LU分解法等算法;然后,对给定的线性方程组进行求解,记录计算结果;最后,分析算法的稳定性和误差。
2. 矩阵特征值问题计算实验:编写程序实现幂法、逆幂法等算法;然后,对给定的矩阵进行特征值和特征向量的计算,记录计算结果;最后,分析算法的收敛性和精度。
3. 非线性方程求根实验:编写程序实现二分法、牛顿法、不动点迭代法等算法;然后,对给定的非线性方程进行求根,记录计算结果;最后,比较不同算法的优缺点。
4. 函数插值实验:编写程序实现拉格朗日插值、牛顿插值等方法;然后,对给定的函数进行插值,记录计算结果;最后,分析插值误差。
5. 常微分方程初值问题数值解法实验:编写程序实现欧拉法、改进的欧拉法、龙格-库塔法等算法;然后,对给定的常微分方程初值问题进行求解,记录计算结果;最后,比较不同算法的稳定性和精度。
清华大学高等数值分析实验设计及答案
高等数值分析实验一工物研13 成彬彬2004310559一.用CG,Lanczos和MINRES方法求解大型稀疏对称正定矩阵Ax=b作实验中,A是利用A= sprandsym(S,[],rc,3)随机生成的一个对称正定阵,S是1043阶的一个稀疏阵A= sprandsym(S,[],0.01,3);检验所生成的矩阵A的特征如下:rank(A-A')=0 %即A=A’,A是对称的;rank(A)=1043 %A满秩cond(A)= 28.5908 %A是一个“好”阵1.CG方法利用CG方法解上面的线性方程组[x,flag,relres,iter,resvec] = pcg(A,b,1e-6,1043);结果如下:Iter=35,表示在35步时已经收敛到接近真实xrelres= norm(b-A*x)/norm(b)= 5.8907e-007为最终相对残差绘出A的特征值分布图和收敛曲线:S=svd(A); %绘制特征值分布subplot(211)plot(S);title('Distribution of A''s singular values');;xlabel('n')ylabel('singular values')subplot(212); %绘制收敛曲线semilogy(0:iter,resvec/norm(b),'-o');title('Convergence curve');xlabel('iteration number');ylabel('relative residual');得到如下图象:为了观察CG方法的收敛速度和A的特征值分布的关系,需要改变A的特征值:(1).研究A的最大最小特征值的变化对收敛速度的影响在A的构造过程中,通过改变A= sprandsym(S,[],rc,3)中的参数rc(1/rc为A的条件数),可以达到改变A的特征值分布的目的:通过改变rc=0.1,0.0001得到如下两幅图以上三种情况下,由收敛定理2.2.2计算得到的至多叠代次数分别为:48,14和486,由于上实验结果可以看出实际叠代次数都比上限值要小较多。
数值分析实验三(解线性方程组)
《数值分析》实验报告实验编号:实验三课题名称:解线性方程组一、算法介绍1、定义四个函数分别为:DieDai(),Newton(),XianWei(),DuiFen()。
在主函数中输入要选用方法所对应的序号后,用switch语句对函数进行调用。
2、迭代法的主要思想为:保留一个变量在等号的左边,其他都移到等号右边。
3、Newton法的主要算法为:把f(x)在x0附近展开成Taylor级数,取其线性部分,选取一点x0,该点所对应的值为f(x0),对于n=1,2,…,Nmax,按Xn+1=Xn-f(Xn)/f’(Xn)求出Xn+1,并计算f(Xn+1),若|Xn+1-Xn|小于容许误差,则停止计算。
4、弦位法:选定初始值x0,x1,并计算f(x0),f(x1);按迭代公式xn+1=xn-f(xn)(xn-xn-1)/(f(xn)-f(xn-1))计算x2,再求f(x2);如果相邻两次迭代值之差在容许的误差范围内则迭代停止,否则用(x2,f(x2)),(x1,f(x1))分别代替(x1,f(x1)),(x0,f(x0)),重复前两个步骤,直至相邻两次迭代值之差在容许的误差范围内。
5、对分区间法:选取方程的根所在的区间(a,b),取其中点c代入方程中得其值为f(c),如果f(c)与f(a)异号说明方程的根在(a,c)区间中,则令b=c,否则令a=c,如果f(c)的绝对值小于0.0001,则停止运算,否则继续计算,直至f(c)的绝对值小于0.0001。
二、程序代码#include <iostream>#include <iomanip>#include <cmath>using namespace std;double f(double x){x=x*x*x-3*x-1;return x;}void DieDai(){cout<<"迭代序列:\n";double x0=2,x1=0;int i=1;while(fabs(x0-x1)>=0.0001){x1=x0;x0=pow((3*x0+1),1.0/3);cout<<"x"<<i<<"="<<setiosflags(ios::fixed)<<setprecision(4)<<x0<<",f(x"<<i<<")="<<f(x0)<<",误差为:"<<x0*x0*x0-3*x0-1<<endl;i++;}cout<<"方程近似解x*="<<x0<<endl;cout<<"共进行"<<i-1<<"次迭代\n" ;}void Newton(){cout<<"迭代序列:\n";double x0=2,x1=0;int i=1;while(fabs(x0-x1)>=0.0001){x1=x0;x0=x0-f(x0)/(3*x0*x0-3);cout<<"x"<<i<<"="<<setiosflags(ios::fixed)<<setprecision(4)<<x0<<",f(x"<<i<<")="<<f(x0)<<",误差为:"<<x0*x0*x0-3*x0-1<<endl;i++;}cout<<"方程近似解x*="<<x0<<endl;cout<<"共进行"<<i-1<<"次迭代\n" ;}void XianWei(){double x0=1,x1=3,x2=2;cout<<"选定曲线y=f(x)上的两个点P0("<<x0<<","<<f(x0)<<")和P1("<<x1<<","<<f(x1)<<")\n";int i=2;while(fabs(f(x2))>=0.0001){x2=x1-f(x1)*(x1-x0)/(f(x1)-f(x0));cout<<setiosflags(ios::fixed)<<setprecision(4)<<"当前区间(";cout<<x0<<","<<x1<<"),与x轴交点("<<x2<<","<<f(x2)<<"),误差为:";if(f(x2)*f(x0)<0){cout<<x2*x2*x2-3*x2-1<<endl;x1=x2;}else{cout<<x2*x2*x2-3*x2-1<<endl;x0=x2;}i++;}cout<<"方程近似解x*="<<x2<<endl;cout<<"共进行"<<i-2<<"次迭代\n" ;}void DuiFen(){double a=1,b=3,c;cout<<"f(x)=0的根的存在区间("<<a<<","<<b<<")\n";cout<<"端点函数值f(a)="<<f(a)<<",f(b)="<<f(b)<<endl;int i=2;while(fabs(f(c))>=0.0001){c=(a+b)/2;cout<<setiosflags(ios::fixed)<<setprecision(4)<<"当前区间("<<a<<","<<b<<"),区间中点x="<<c;cout<<",f(x)="<<f(c)<<",误差为:";if(f(c)*f(a)<0){cout<<c*c*c-3*c-1<<endl;b=c;}else{cout<<c*c*c-3*c-1<<endl;a=c;}i++;}cout<<"方程近似解x*="<<c<<endl;cout<<"共进行"<<i-2<<"次迭代\n" ;}void Menu(){int n;cin>>n;switch(n){case 1:cout<<"*** 迭代法***\n";DieDai();Menu();break;case 2:cout<<"*** Newton法***\n";Newton();Menu();break;case 3:cout<<"*** 弦位法***\n";XianWei();Menu();break;case 4:cout<<"*** 对分区间法***\n";DuiFen();Menu();break;case 5:return;}}int main (){cout<<" *****问题:求f(x)=x*x*x-3x-1=0在x0=2附近的实根。
【清华】7.0_实验5-线性代数方程组的数值解法
0.99995
1
1
0
0.5833
0.8333
0.9967
0.99992
1
1
0
0.5000
0.7847
0.9954
0.99988 0.99999
1
0
0.5000
0.7569
0.9943
0.99985 0.99999
1
0
0.5000
0.7500
0.9935
0.99983 0.99999
1
0
0.5000
次时得到的向量,可以作出论断,G-S 方法已经基本可以认为算到了精确值 (1,1,,1)T (共
20 个 1),而此时 JACOBI 法也接近精确值。由此说明,高斯-赛德尔方法比雅可比方法具有 更快的收敛速度。原因就在于 G-S 法每步都利用了当前步骤计算出的最新的 x 分量值。
进一步,再选取不同的 x0 和 b 来验证:
0.7500
0.9930
0.99981 0.99999 0.99999
0
0.5000
0.7500
0.9926
0.99979 0.99999 0.99999
0
0.5000
0.7500
0.9924
0.99978 0.99999 0.99999
0
0.5000
0.7500
0.9923
0.99977 0.99999 0.99999
序如下:
b0=ones(20,1);x0=zeros(20,1);
%给 b 和 x0 赋值
a1=A;
x11=jacobi(x0,b0,a1,m); n1=size(x11)
大学数学实验五_线性代数方程组的数值解法
【实验目的】 1、学会用 MATLAB 软件数值求解线性代数方程组,对迭代法的收敛性和解
的稳定性作初步分析。 2、通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】
3 已知方程组 Ax=b,其中
,定义为
试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对 收敛速度的影响。实验要求: (1) 选取不同的初始向量 x(0)和不同的方程组的右端项向量 b,给定迭代误差要求,用雅
k=k+1; xj=Bj*xj+fj; 多输出了矩阵 P,矩阵 P 可视为一个行向量,其每个元素均为迭代 k 次后得到的 xk。这样以 k 为横轴,解向量为纵轴,可输出图形观察 xk 是否收敛。函数 GaussSeidel 也需作同样修改,修改后的函数在此不再赘述。
模型: 已知某年该植物的数量为 x0,记第 k 年的植物数量为 xk,那么有 xk + pxk-1 + qxk-2 = 0 (k = 2, 3, …… , n)
其中 p = -a1bc,q = -a2b(1-a1)bc。若要求 n 年后数量达到 xn,则 Ax = b
其中
,
,
7
① 用稀疏系数矩阵求解。
这个函数中,n 表示矩阵 A 的阶数,在本题中恒取 20,a 表示主对角线元素的值,b 在 本题中恒取-1/4,c 在本题中恒取-1/2。
编写用雅可比迭代法求方程解的函数 Jacobi。
function [xj,k]=Jacobi(A,X0,b,e) D=diag(diag(A)); n=length(A); L=-(tril(A)-D); U=-(triu(A)-D); fj=D\b; Bj=D\(L+U); xj=X0; k=0; while norm(A*xj-b)/norm(b)>e
【数值分析实验报告感想】数值分析实验报告
【数值分析实验报告感想】数值分析实验报告数值分析实验报告数值分析实验报告学院:姓名:实验一一、实验名称:线性方程组的直接解法二、实验目的:1、理解Gauss列主元素法和LU分解法的基本原理2、掌握Gauss列主元素法和LU分解法3、编写maTLaB程序实现Gauss列主元素法和LU分解法三、实验步骤及内容:使用Gauss列主元素法和LU分解法求解下列方程组:5.0x1?0.96x2?6.5x3?0.962.0x1?4.5x2?0.36x3?0.0200.50x1?1.1x2?3.1x3?6.01.基本计算公式:(一)Gauss列主元素法:主要包括消元和回代步骤:[1]消元过程:主要是把原方程组化为三角方阵的方程组也即lik?aik/akkaijk?1kkkbik?1?bik?likbk?aij?likakj与Gauss消元法的区别在于在每次消元前寻找该列绝对值最大的元素作为主元,然后进行消元。
[2]回代过程:按变量的逆序逐步回代可得到方程组的解:有xn?b/axk?(b?nnnnnkkkkax)/a?kllkk(n-1,n-2,?2,1,)nl?k?1根据以上公式可以求出方程的解。
(二)LU分解法:其基本思想是:对系数矩阵a进行分解,a?L?U,分别得到L为下三角矩阵,U为上三角矩阵,最后通过两个三角方程组Ly?b和Ux?y,最后得出方程组的解。
2、基本程序设计:列主元素法消元法的maTLaB基本程序如下所示:function[x,det,flag]=Gauss(a,b)%求线性方程组的列主元Gauss消元法,其中%a为方程组的系数矩阵;%b为方程组的右端项;%x为方程组的解;%det为系数矩阵a的行列式的值;%flag为指标向量,flag='failure'表示计算失败,flag=‘oK’表示计算成功。
a=[5.00.966.5;2.04.50.36;0.501.13.1];b=[0.960.0206.0];[n,m]=size(a);nb=length(b)%当方程组行与列的维数不相等时,停止计算,并输出出错信息ifn~=merror('Therowandcolumnsofmatrixamustbeepual!'); return;end%当方程组与右端项的维数不匹配时,停止计算,并输出错误信息ifn~=nberror('Thecolumnsofamustbeequalthelengthofb!')return;end%开始计算,先赋初值flag='oK';det=1;x=zeros(n,1);fork=1:n-1%选主元max1=0;fori=k:nifabs(a(i,k))>max1max1=abs(a(i,k));r=i;endendifmax1flag='failure';return; end%交换两行ifr>kforj=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;det=-det;end%消元过程fori=k+1:nm=a(i,k)/a(k,k);forj=k+1:na(i,j)=a(i,j)-m*a(k,j);endb(i)=b(i)-m*b(k);enddet=det*a(k,k);enddet=det*a(n,n);%回代过程ifabs(a(n,n))flag='failure';return; endfork=n:-1:1forj=k+1:nb(k)=b(k)-a(k,j)*x(j);endx(k)=b(k)/a(k,k);endx(k)=b(k)/a(k,k);endLU分解法求解线性方程组程序代码:functiony=isexist(a,b)[m,n]=size(a);[m1,n1]=size(b);ifm~=m1error('您的输入有误!!!'); return;endr=rank(a);s=rank([a,b]);ifr==s&r==ny=1;elseifr==s&ry=inf;elsey=0;endfunctionx=LUSeparation(a,b)flag=isexist(a,b);ifflag==0disp('该线性方程组无解!!!'); x==[];return;elser=rank(a);[m,n]=size(a);[L,U,p]=lu(a);b=p*b;y(1)=b(1);ifm>1fori=2:my(i)=b(i)-L(i,1:i-1)*y(1:i-1)';endendy=y';x0(r)=y(r)/U(r,r);ifr>1fori=r-1:-1:1x0(i)=(y(i)-U(i,i+1:r)*x0(i+1:r)')/U(i,i); endendx0=x0';ifflag==1x=x0;return;elseformatrat;z=null(a,'r');[m2,n2]=size(z);x0(r+1:n)=0;fori=1:n2;t=sym(char([10748+i]));k(i)=t;endx=x0;fori=1:n2x=x+k(i)*z(:,i);endendend下页数值分析实验报告四、实验结果:由列主元素消元法解的方程组的解为:在命令窗口中输入:LU分解法运行结果:五、总结:列主元素法相比高斯消元法而言,其避免了小主元的出现的问题,抑制了舍入误差的增长,使得计算结果相比而言误差较小。
【清华】实验5-线性代数方程组的数值解法(011813)
实验5 线性代数方程组的数值解法化学工程系分9班焦阳2009011813 【实验目的】1. 学会用MATLAB软件数值求解线性代数方程组,对迭代法的收敛性和解的稳定性作初步分析;2. 通过实例学习用线性代数方程组解决简化的实际问题。
【实验内容】1.题目3:已知方程组,其中,定义为:试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方程组系数矩阵性质对收敛速度的影响。
实验要求:(1)选取不同的初始向量和不同的方程组右端向量b,给出迭代误差要求,用雅可比迭代法和高斯-赛德尔迭代法计算,观测得到的迭代向量序列是否均收敛?若收敛,记录迭代次数,分析计算结果并得出你的结论;(2)取定右端向量b和初始向量,将A的主对角线元素成倍增长若干次,非主对角线元素不变,每次用雅可比迭代法计算,要求迭代误差满足 ,比较收敛速度,分析现象并得出你的结论。
迭代法的原理分析:对于线性方程组,满足A的对角元素均不为0,则可以将A分解为:。
则雅可比迭代公式即为:,又由于:因此有:,同理,高斯-赛德尔迭代法的矩阵表示为:,根据原理,先编写雅可比迭代法和高斯-赛德尔迭代法的函数程序。
(1) 一、雅可比迭代: I.选取初始向量,。
迭代精度的要求是,运用雅可比迭代进行求解,每次迭代后的结果如下表:迭代次数0 1 23 456 x1 1 0.583333 0.527778 0.499421 0.489583 0.485139 0.483239 x2 1 0.75 0.638889 0.602431 0.586034 0.579162 0.576046 x3 1 0.833333 0.715278 0.670139 0.649595 0.640545 0.636396 x4 1 0.833333 0.743056 0.693287 0.671586 0.661185 0.6564 x5 1 0.833333 0.75 0.704282 0.681713 0.670915 0.665718 x6 1 0.833333 0.75 0.707755 0.685571 0.674716 0.669367 x7 1 0.833333 0.75 0.708333 0.687066 0.676272 0.670924 x8 1 0.833333 0.75 0.708333 0.687452 0.67685 0.671535 x9 1 0.833333 0.75 0.708333 0.6875 0.677039 0.671767 x10 1 0.833333 0.75 0.708333 0.6875 0.677079 0.671844 x11 1 0.833333 0.75 0.708333 0.6875 0.677079 0.671844 x12 1 0.833333 0.75 0.708333 0.6875 0.677039 0.671767 x13 1 0.833333 0.75 0.708333 0.687452 0.67685 0.671535 x14 1 0.833333 0.75 0.708333 0.687066 0.676272 0.670924 x15 1 0.833333 0.75 0.707755 0.685571 0.674716 0.669367 x16 1 0.833333 0.75 0.704282 0.681713 0.670915 0.665718 x17 1 0.833333 0.743056 0.693287 0.671586 0.661185 0.6564 x18 1 0.833333 0.715278 0.670139 0.649595 0.640545 0.636396 x19 1 0.75 0.638889 0.602431 0.586034 0.579162 0.576046 x20 1 0.5833330.5277780.499421 0.489583 0.485139 0.483239 迭代次数78 9 10 11 12 x1 0.48237396 0.48198 0.481798 0.481713 0.481672 0.481653 x2 0.57463907 0.573988 0.573686 0.573543 0.573476 0.573444 x3 0.63448738 0.633597 0.633179 0.632981 0.632887 0.632843 x40.65413674 0.653071 0.6525650.6523240.6522080.652153x5 0.66323779 0.662048 0.661477 0.661202 0.661069 0.661005x6 0.66676822 0.665504 0.66489 0.664591 0.664446 0.664376x7 0.66827415 0.666972 0.666333 0.666019 0.665865 0.66579x8 0.66888283 0.667565 0.666913 0.66659 0.666431 0.666353x9 0.66912708 0.667806 0.667148 0.666821 0.666659 0.666578x10 0.66921041 0.66789 0.66723 0.666901 0.666738 0.666657x11 0.66921041 0.66789 0.66723 0.666901 0.666738 0.666657x12 0.66912708 0.667806 0.667148 0.666821 0.666659 0.666578x13 0.66888283 0.667565 0.666913 0.66659 0.666431 0.666353x14 0.66827415 0.666972 0.666333 0.666019 0.665865 0.66579x15 0.66676822 0.665504 0.66489 0.664591 0.664446 0.664376x16 0.66323779 0.662048 0.661477 0.661202 0.661069 0.661005x17 0.65413674 0.653071 0.652565 0.652324 0.652208 0.652153x18 0.63448738 0.633597 0.633179 0.632981 0.632887 0.632843x19 0.57463907 0.573988 0.573686 0.573543 0.573476 0.573444x20 0.48237396 0.48198 0.481798 0.481713 0.481672 0.481653 通过表1可知,一共迭代了12次。
数学实验 5:线性代数方程组的数值解法
大倍数 1.0000 2.0000 3.0000 4.0000 5.0000 6.0000 7.0000 8.0000 9.0000 10.0000 11.0000 12.0000 13.0000 14.0000 15.0000 16.0000 17.0000 18.0000 19.0000 20.0000
小值q 0.4893 0.2447 0.1631 0.1223 0.0979 0.0816 0.0699 0.0612 0.0544 0.0489 0.0445 0.0408 0.0376 0.0350 0.0326 0.0306 0.0288 0.0272 0.0258 0.0245 21.0000 12.0000 9.0000 8.0000 8.0000 7.0000 7.0000 7.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 6.0000 5.0000 5.0000
实验 5:线性代数方程组的数值解法
习题3:
已知方程组,其中,定义为: 试通过迭代法求解此方程组,认识迭代法收敛的含义以及迭代初值和方 程组系数矩阵性质对收敛速度的影响。实验要求: (1) 选取不同的初始向量x0和不同的方程组右端向量b,给定迭 代误差要求,用雅可比迭代法和高斯-赛德尔迭代法计算, 观测得到的迭代向量序列是否均收敛?若收敛,记录迭代 次数,分析计算结果并得出结论; (2) 取定右端向量b和初始向量x0,将A的主对角线元素成倍的 增长若干次,非主对角元素不变,每次用雅可比迭代法计 算,要求迭代误差满足,比较收敛速度,分析现象并得出结 论。 1、 程序设计(可直接粘贴运行) 1) Jacobi迭代法 function y=jacobi(a,b,x0,e,m) %定义jacobi函数,其中:a,b为线性方程组中的矩阵和右端向量;x0 为初始值; %e和m分别为人为设定的精度和预计迭代次数;运行结果y为迭代的结 果和所有中间值组成的 %矩阵 y=0; %对y初始化 d=diag(diag(a)); %按雅可比迭代标准形 形式取主对角元素作为矩阵D u=-triu(a,1); %取上三角矩阵u l=-tril(a,-1); %取下三角矩阵l bj=d^-1*(l+u); fj=d^-1*b; x=[x0,zeros(20,m-1)]; %初始化x,其中x1=x0,即 初始值 for k=1:m %人为规定迭代次 数,防止不收敛迭代导致死循环 x(:,k+1)=bj*x(:,k)+fj; %jacobi迭代 if norm(x(:,k+1)-x(:,k),inf)<e
数值分析的实验报告
数值分析的实验报告数值分析的实验报告导言数值分析是一门研究数值计算方法和数值计算误差的学科,它在科学计算、工程技术和社会经济等领域具有广泛的应用。
本实验旨在通过对数值分析方法的实际应用,验证其有效性和可靠性。
实验一:方程求根方程求根是数值分析中的基础问题之一。
我们选取了一个非线性方程进行求解。
首先,我们使用二分法进行求解。
通过多次迭代,我们得到了方程的一个近似解。
然后,我们使用牛顿法进行求解。
与二分法相比,牛顿法的收敛速度更快,但需要选择一个初始点。
通过比较两种方法的结果,我们验证了牛顿法的高效性。
实验二:插值与拟合插值与拟合是数值分析中常用的数据处理方法。
我们选取了一组实验数据,通过拉格朗日插值法和最小二乘法进行插值和拟合。
通过对比两种方法的拟合效果,我们验证了最小二乘法在处理含有噪声数据时的优势。
同时,我们还讨论了插值和拟合的精度与样本点数量之间的关系。
实验三:数值积分数值积分是数值分析中的重要内容之一。
我们选取了一个定积分进行计算。
首先,我们使用复化梯形公式进行积分计算。
通过增加分割区间的数量,我们得到了更精确的结果。
然后,我们使用复化辛普森公式进行积分计算。
与复化梯形公式相比,复化辛普森公式具有更高的精度。
通过比较两种方法的结果,我们验证了复化辛普森公式的优越性。
实验四:常微分方程数值解常微分方程数值解是数值分析中的重要应用之一。
我们选取了一个常微分方程进行数值解的计算。
首先,我们使用欧拉方法进行数值解的计算。
然后,我们使用改进的欧拉方法进行数值解的计算。
通过比较两种方法的结果,我们验证了改进的欧拉方法的更高精度和更好的稳定性。
实验五:线性方程组的数值解法线性方程组的数值解法是数值分析中的重要内容之一。
我们选取了一个线性方程组进行数值解的计算。
首先,我们使用高斯消元法进行数值解的计算。
然后,我们使用追赶法进行数值解的计算。
通过比较两种方法的结果,我们验证了追赶法在求解三对角线性方程组时的高效性。
工程数值分析实验报告(3篇)
第1篇一、实验目的本次实验旨在通过数值分析的方法,对工程实际问题进行建模、求解和分析。
通过学习数值方法的基本原理和算法,提高解决实际工程问题的能力。
二、实验内容1. 线性方程组的求解2. 矩阵特征值与特征向量的计算3. 函数插值与曲线拟合4. 数值微分与积分三、实验步骤1. 线性方程组的求解(1)编写程序实现高斯消元法、克劳斯消元法和列主元素法(2)设计输入界面,用户输入增广矩阵的行和列,填写系数及常数项(3)分别运用三种方法求解线性方程组,比较求解结果的正确性、数值稳定性和计算效率2. 矩阵特征值与特征向量的计算(1)编写程序实现幂法、QR算法和逆幂法(2)设计输入界面,用户输入矩阵的行和列,填写矩阵元素(3)分别运用三种方法计算矩阵的特征值与特征向量,比较求解结果的准确性和计算效率3. 函数插值与曲线拟合(1)编写程序实现拉格朗日插值、牛顿插值和样条插值(2)设计输入界面,用户输入函数的自变量和函数值,选择插值方法(3)分别运用三种方法进行函数插值,比较插值结果的准确性和光滑性4. 数值微分与积分(1)编写程序实现有限差分法、龙格-库塔法和辛普森法(2)设计输入界面,用户输入函数的导数或积分的上下限,选择数值方法(3)分别运用三种方法进行数值微分和积分,比较求解结果的准确性和计算效率四、实验结果与分析1. 线性方程组的求解通过实验,我们发现列主元素法在求解线性方程组时具有较好的数值稳定性,计算效率也较高。
而高斯消元法和克劳斯消元法在处理大型稀疏矩阵时存在一定的困难。
2. 矩阵特征值与特征向量的计算实验结果表明,QR算法和逆幂法在计算矩阵特征值与特征向量时具有较高的准确性和计算效率。
幂法在处理大型稀疏矩阵时表现出较好的性能。
3. 函数插值与曲线拟合在函数插值和曲线拟合实验中,样条插值方法具有较好的准确性和光滑性。
拉格朗日插值和牛顿插值方法在处理简单函数时表现良好,但在处理复杂函数时可能存在精度问题。
数值分析实验报告
实验五 解线性方程组的直接方法实验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)扰动相对误差与解的相对偏差,分析它们与条件数的关系。
数值分析实验报告_清华大学__线性代数方程组的数值解法
线性代数方程组的数值解法实验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)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
1.1程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A); Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx1.2运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =2.5575×103Cond(A,2) = 1.7276×103Cond(A,inf) =2.5575×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T (2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3 …(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T (3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000,1.000000000000000,1.000000000000000,1.000000000000001,0.999999999999998,1.000000000000004,0.999999999999993,1.000000000000014,0.999999999999972,1.000000000000057,0.999999999999886,1.000000000000227,0.999999999999547,1.000000000000902,0.999999999998209,1.000000000003524,0.999999999993179,1.000000000012732,0.999999999978173,1.000000000029102]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵将计算结果列于下表:1.3简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
清华大学数值分析课件解线性代数方程组的直接法(2)
l13 l23
l1n l2 n 1
ln1 1
tik lik d k
ki
14
改进平方根法(无开方Choleskey分解)解对称 正定方程组的算法
d1 a11 对于 i 2, 3, , n, 计算 j 1 lij aij lik d k l jk d j , k 1 i 1 2 di aii d k lik k 1
20
第三步:由 L y f 先求 y: 1 l2 Ly f y1 y 2 y 3 y4 1 l3 y1 y y2 3 y 1 4 y1 f1 f1 l y y f 2 2 2 1 f2 f3 l 3 y2 y3 f 3 f 4 l y y f 4 4 4 3
一般地:
对 i 3,
ri ri
, n 计算
r 1 k 1 rk ki
u a l u l
ir
ir ,
,n ,n
(a ir
k 1
r 1
l u )/u
ik kr
rr
i r 1 ,
注:该公式的特点:U 的元素按行求, L 的元素按列求; 先求
U 的第k 行, 再求 L 的第k 列, U 和 L 一行一列交叉计算.
10
用平方根法解线性方程组Ax = b的算法
(1)对矩阵A进行Cholesky分解,即A=LLT,由矩阵乘法: 对于 j= 1, 2,…, n 计算
j 1 2 l jj a jj l jk k 1 1 2
j 1 lij aij lik l jk l jj k 1
实验6线性代数方程组的数值解法
实验6线性代数⽅程组的数值解法实验6 线性代数⽅程组的数值解法[实验⽬的]1. 1. 学会⽤MATLAB 软件数值求解线性代数⽅程组,对迭代法的收敛性和解的稳定性作初步分析;2. 2. 通过实例学习⽤线性代数⽅程组解决简化的实际问题。
[实验内容]5-5 输电⽹络:⼀种⼤型输电⽹络可简化为图5.5(见书)所⽰电路,其中R 1,R 2,…,R n 表⽰负载电阻,r 1,r 2,…,r n 表⽰线路内阻,I 1,I 2,…,I n 表⽰负载上的电流。
设电源电压为V 。
(1)列出求各负载电阻R 1,R 2,…,R n 的⽅程;(2)设I 1=I 2=…=I n =I ,r 1=r 2=…=r n =r ,在r=1,I=0.5,V=18,n=10的情况下求R 1,R 2,…,R n 及总电阻R 0。
[问题分析、模型建⽴及求解](1) 设电源负极为电势为0,电阻R 1上对应节点电压为V 1,对于任意节点,根据KCL 定律列出⽅程:111++----=k k k k k k k k r V V r V V R V⽽k kk R V I =,可得:111111)(++++--++-=k k k k k k k k k k k k R r IR r I r I R r I Ik=2,3,……,n-1;k=1时2221211R r I R r I I +-=,为与上式形式⼀致,化为22212111111)(R r IR r I r I r V I +--=-k=m (12-≤≤n m )时 111111)(++++--+--+=m m m n m m m m m m m m R r IR r I r I R r I Ik=n 时R r I I -=--11设以上⽅程组的矩阵形式为:b AR = 则 []Tn R R R R ΛΛΛΛ21=Tn I I I r V I b ?-=ΛΛΛ3211---------=------n n nn n n nn n n n n r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I r I A 11 11114443333233322221222111000000ΛΛΛO O M O O O O M M O O O O O M M O M O ΛΛ5.021=====I I I I n Λ,121=====r r r r n Λ,V=18,n=10,-----=5.05.0005.015.005.015.005.015.0005.01ΛΛΛO M O O O O M M O O O O O MM O M O ΛΛΛA[]Tb 5.05.05.05.17ΛΛΛ-=在命令窗⼝输⼊MA TLAB 程序如下:clear all ;n=10; %由题⽬要求设定A11=sparse(1:n-1,1:n-1,-1,n,n); %定义A 的对⾓元素,除(n,n) A12=sparse(n,n,-0.5,n,n); %定义(n,n)A1=A11+A12; %对⾓元素A2=sparse(1:n-1,2:n,0.5,n,n); %输⼊A 的上次对⾓元素 A3=sparse(2:n,1:n-1,0.5,n,n); %输⼊A 的下次对⾓元素 A=A1+A2+A3; b1=0.5*ones(n,1); %b 的除第⼀项元素 b2=sparse(1,1,18,n,1); %b 的第⼀项元素 b=b1-b2; R=A\b输出结果如下:R =26.000017.0000 9.0000 2.0000 -4.0000 -9.0000 -13.0000 -16.0000 -18.0000 -19.0000所以各阻值为(R 1,R 2,…,R 10)=(26,17,9,2,-4,-9,-13,-16,-18,-19)总电阻R 0(即输⼊等效电阻)为00I V R =,⼜nI得到 )(6.35.010180Ω=?=R5-6 有5个反应器连接如图5.6(见书),各个Q 表⽰外部输⼊、输出及反应器间的流量(m 3/min ),各个c 表⽰外部输⼊及反应器内某物质的浓度(mg/m 3)。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
线性代数方程组的数值解法实验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)选取其他你感兴趣的问题或者随机生成矩阵,计算其条件数。
重复上述实验,观察记录并分析实验结果。
1.1程序清单n=input('矩阵A 的阶数:n=');A=6*diag(ones(1,n))+diag(ones(1,n-1),1)+8*diag(ones(1,n-1),-1); b=A*ones(n,1);p=input('计算条件数使用p-范数,p='); cond_A=cond(A,p) [m,n]=size(A); Ab=[A b];r=input('选主元方式(0:自动;1:手动),r=');Abfor i=1:n-1switch rcase(0)[aii,ip]=max(abs(Ab(i:n,i)));ip=ip+i-1;case (1)ip=input(['第',num2str(i),'步消元,请输入第',num2str(i),'列所选元素所处的行数:']);end;Ab([i ip],:)=Ab([ip i],:);aii=Ab(i,i);for k=i+1:nAb(k,i:n+1)=Ab(k,i:n+1)-(Ab(k,i)/aii)*Ab(i,i:n+1);end;if r==1Abendend;x=zeros(n,1);x(n)=Ab(n,n+1)/Ab(n,n);for i=n-1:-1:1x(i)=(Ab(i,n+1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endx1.2运行结果(1)n=10,矩阵的条件数及自动选主元Cond(A,1) =2.5575×103Cond(A,2) = 1.7276×103Cond(A,inf) =2.5575×103程序自动选择主元(列主元)a.输入数据矩阵A的阶数:n=10计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=0b.计算结果x=[1,1,1,1,1,1,1,1,1,1]T (2)n=10,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000, 1.000000000000000, 1.000000000000000, 1.000000000000001, 0.999999999999998, 1.000000000000004, 0.999999999999993, 1.000000000000012, 0.999999999999979, 1.000000000000028]Tb. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=10 计算条件数使用p-范数,p=1选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3 …(实际选择时,第k 步选择主元处于第k+1行) 最终计算得x=[1,1,1,1,1,1,1,1,1,1]T (3)n=20,手动选主元a. 每步消去过程总选取按模最小或按模尽可能小的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:1(2)(2) 6.0000 1.00007.00004.6667 1.0000 5.66678.0000 6.000015.0000[]8.00001.000015.00006.0000 1.00008.0000 6.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:2 …(实际选择时,第k 步选择主元处于第k 行) 最终计算得x=[1.000000000000000,1.000000000000000,1.000000000000000,1.000000000000001,0.999999999999998,1.000000000000004,0.999999999999993,1.000000000000014,0.999999999999972,1.000000000000057,0.999999999999886,1.000000000000227,0.999999999999547,1.000000000000902,0.999999999998209,1.000000000003524,0.999999999993179,1.000000000012732,0.999999999978173,1.000000000029102]T b. 每步消去过程总选取按模最大的元素作为主元 矩阵A 的阶数:n=20 计算条件数使用p-范数,p=1 选主元方式(0:自动;1:手动),r=1(1)(1)61786115[]861158614A b ⎡⎤⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦第1步消元,请输入第1列所选元素所处的行数:2(2)(2)8.0000 6.0000 1.000015.0000-3.50000.7500-4.250008.0000 6.0000 1.000015.0000[]8.0000 6.000015.00008.0000 1.00006.0000 1.000015.00008.0000 6.000014.0000A b ⎡⎤⎢⎥-⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第2步消元,请输入第2列所选元素所处的行数:3…(实际选择时,第k步选择主元处于第k+1行)最终计算得x=[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]T(4)A分别为幻方矩阵,Hilbert矩阵,pascal矩阵和随机矩阵将计算结果列于下表:1.3简要分析计算(1)表明:对于同一矩阵,不同范数定义的条件数是不同的;Gauss消去法在消去过程中选择模最大的主元能够得到比较精确的解。
计算(2)表明:通过比较每次选取模最大或模最小的元素的选主元方式,可以发现,在本题给定的问题中,选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确。
因为这样做可以避免使用较小的数做为除数,以免发生结果数量级加大,使大数吃掉小数,产生舍入误差。
计算(3)表明:首先,n=20得到与计算(2),即n=10时一样的结论,即选取模最大的元素作为主元比取模最小的元素作为主元时产生的结果更精确;其次,与计算(2) (Cond(A10×10,1) =2.5575×103)比较,Cond(A20×20,1) =2.6214×106显著增大,且计算(3)的误差也远大于计算(2),即,矩阵的条件数越大,产生的误差也越大。
计算(4)表明: Gauss消去法在消去过程中,主元的选择与算法的稳定性有密切的联系。
一般来说,选取绝对值大的元素作为主元比绝对值小的元素作为主元时的计算结果更加精确。
但这并不是绝对的,一些特殊的方阵,如Pascal方矩阵,则恰恰是选择模最小的元素作为主元时计算结果最精确(选模最小的元素只是一个表象,这种选主元方法优于其他选主元方法的本质是这种选择方法能使消去过程不产生浮点数,而全是整数运算,只有在回代过程中才有可能会产生浮点数)。
在系数矩阵性质未知,或者说对于绝大多数的系数矩阵来说,选择模最大的元素作为主元是一种比较稳定和精确的方法。
实验2.病态的线性方程组的求解问题提出:理论的分析表明,求解病态的线性方程组是困难的。
实际情况是否如此,会出现怎样的现象呢?实验内容:考虑方程组Hx=b 的求解,其中系数矩阵H 为Hilbert 矩阵,n j i j i h h H j i n n j i ,,2,1,,11,)(,, =-+==⨯这是一个著名的病态问题。