MATLAB数值分析实验三(线性方程求解及精度分析)

合集下载

数值分析中求解线性方程组的MATLAB程序(6种)

数值分析中求解线性方程组的MATLAB程序(6种)

数值分析中求解线性方程组的MATLAB程序(6种)1.回溯法(系数矩阵为上三角)function X=uptrbk(A,B)%求解方程组,首先化为上三角,再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));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;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵(先化为上三角,在用回溯法)function X=uptrbk(A,B)%求解方程组,首先化为上三角,再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));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;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B,L为下三角,U为上三角%程序中并没有真正解出L和U,全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for 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=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=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);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法,求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。

MATLAB计算方法3解线性方程组计算解法

MATLAB计算方法3解线性方程组计算解法

MATLAB计算方法3解线性方程组计算解法线性方程组是数学中的一个重要问题,解线性方程组是计算数学中的一个基本计算,有着广泛的应用。

MATLAB是一种功能强大的数学软件,提供了多种解线性方程组的计算方法。

本文将介绍MATLAB中的三种解线性方程组的计算方法。

第一种方法是用MATLAB函数“linsolve”解线性方程组。

该函数使用高斯消元法和LU分解法求解线性方程组,可以处理单个方程组以及多个方程组的情况。

使用该函数的语法如下:X = linsolve(A, B)其中A是系数矩阵,B是常数向量,X是解向量。

该函数会根据A的形式自动选择求解方法,返回解向量X。

下面是一个使用“linsolve”函数解线性方程组的例子:A=[12;34];B=[5;6];X = linsolve(A, B);上述代码中,A是一个2×2的系数矩阵,B是一个2×1的常数向量,X是一个2×1的解向量。

运行代码后,X的值为[-4.0000;4.5000]。

第二种方法是用MATLAB函数“inv”求解逆矩阵来解线性方程组。

当系数矩阵A非奇异(可逆)时,可以使用逆矩阵求解线性方程组。

使用“inv”函数的语法如下:X = inv(A) * B其中A是系数矩阵,B是常数向量,X是解向量。

该方法先计算A的逆矩阵,然后将逆矩阵与B相乘得到解向量X。

下面是一个使用“inv”函数解线性方程组的例子:A=[12;34];B=[5;6];X = inv(A) * B;上述代码中,A是一个2×2的系数矩阵,B是一个2×1的常数向量,X是一个2×1的解向量。

运行代码后,X的值为[-4.0000;4.5000]。

第三种方法是用MATLAB函数“mldivide”(或“\”)求解线性方程组。

该函数使用最小二乘法求解非方阵的线性方程组。

使用“mldivide”函数的语法如下:X=A\B其中A是系数矩阵,B是常数向量,X是解向量。

数值分析 第三章 基于MATLAB的科学计算—线性方程组1概述

数值分析 第三章 基于MATLAB的科学计算—线性方程组1概述

科学计算—理论、方法及其基于MATLAB的程序实现与分析 三、 解线性方程组(线性矩阵方程)解线性方程组是科学计算中最常见的问题。

所说的“最常见”有两方面的含义:1) 问题的本身是求解线性方程组;2) 许多问题的求解需要或归结为线性方程组的求解。

关于线性方程组B A x B Ax 1-=⇒=(1)其求解方法有两类:1) 直接法:高斯消去法(Gaussian Elimination ); 2) 间接法:各种迭代法(Iteration )。

1、高斯消去法1) 引例考虑如下(梯形)线性方程组:()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路:把一般的线性方程组(1)化成(上或下)梯形的形式。

2)高斯消去法——示例考虑如下线性方程组:⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x1) 第一个方程的两端乘12加到第二个方程的两端,第一个方程的两端乘-1加到第三个方程的两端,得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x2) 第二个方程的两端乘001.010-加到第三个方程的两端,得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x3) 从上述方程组的第三个方程依此求解,得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x 3)高斯消去法的不足及其改进——高斯(全、列)主元素消去法在上例中,由于建模、计算等原因,系数2.001而产生0.0005的误差,实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x 注:数值稳定的算法高斯列主元素消去法就是在消元的每一步选取(列)主元素—一列中绝对值最大的元取做主元素,高斯列主元素消去法是数值稳定的方法。

matlab线性方程组数值求解实验报告

matlab线性方程组数值求解实验报告

湖南大学电气与信息工程学院 《数值计算》课程 上机实验报告一. 实验目的:了解gauss 消去法和迭代法matlab 算法实现求任意方程组的根。

二. 实验内容:用gauss 消去法和迭代法求解下列线性方程组:263234323923321321321=++=++=++x x x x x x x x x1.求出gauss 消去法的上三角矩阵和方程组的解321,,x x x ,并在命令窗口显示;2.显示迭代法求解过程中所有结果(,,,,,,,,,321131*********NN N x x x x x x x x x ⋯⋯)要求求解精度达到10^-5.三. 算法介绍或方法基础1) 消去法:消元过程:设0)0(11≠a ,令乘数)0(11)0(11/a a m i i -=,做(消去第i 个方程组的i x )操作1i m ×第1个方程+第i 个方程(i=2,3,.....n )则第i 个方程变为1)1(2)1(2...i n in i b x a x a =++ 这样消去第2,3,。

,n 个方程的变元i x 后。

原线性方程组变为:⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=++=++)1()1(2)1(2)1(2)1(22)1(22)0(1)0(11)0(11... . .... ...n n nn n n n n n b x a x a b x a x a b x a x a 这样就完成了第1步消元。

回代过程:在最后的一方程中解出n x ,得:)1()1(/--=n nn n n n a b x再将n x 的值代入倒数第二个方程,解出1-n x ,依次往上反推,即可求出方程组的解: 其通项为3,...1-n 2,-n k /)()1(1)1()1(=-=-+=--∑k kk nk j j k kj k kk a x abx高斯赛德尔迭代法:由雅可比迭代公式可知,在迭代的每一步计算过程中是用()k x的全部分量来计算()1+k x的所有分量,显然在计算第i 个分量()1+k ix 时,已经计算出的最新分量()()1111+-+k i k x ,...,x 没有被利用,从直观上看,最新计算出的分量可能比旧的分量要好些.因此,对这些最新计算出来的第1+k 次近似()1+k x 的分量()1+k jx 加以利用,就得到所谓解方程组的高斯—塞德(Gauss-Seidel )迭代法.把矩阵A 分解成U L D A --= (6)其中()nn a ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 ()b Ux x L D +=-即 22f x B x +=其中()()b L D f ,U L D B 1212---=-= (7)以2B 为迭代矩阵构成的迭代法(公式)()()221f x B x k k +=+ (8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++ (9)由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出()1+k ix 后()k ix 不再使用,所以用()1+k ix 冲掉()k ix ,以便存放近似解.四.程序1)消去法:function x=gauss(A,b)n=length(b);A=[A,b];for k=1:(n-1)A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))... -A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1));A((k+1):n,k)=zeros(n-k,1);Aendx=zeros(n,1);x(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1x(k,:)=(A(k,n+1)-A(k,(k+1):n)*x((k+1):n))/A(k,k);end2)迭代法:function EX()a=input('请输入系数矩阵a:');b=input('请输入矩阵b:');N=input('请输入最大迭代次数N:');esp=input('请输入近似解的误差限:');if any(diag(a))==0error('系数矩阵错误,迭代终止!')endD=diag(diag(a));X0=zeros(size(b));x1=0;x2=0;x3=0;X1=[x1;x2;x3];h=inv(D)*b;B=inv(D)*(D-a);B1=triu(B);B2=tril(B);k=1;fprintf('高斯-赛德尔迭代法');fprintf('第0次迭代得:')disp(X1');while k<=Nx1=h(1,1)+B1(1,:)*X0;X1=[x1;x2;x3];x2=h(2,1)+B1(2,:)*X0+B2(2,:)*X1;X1=[x1;x2;x3];x3=h(3,1)+B2(3,:)*X1;X1=[x1;x2;x3];if norm(X1-X0,inf)<espfprintf('已满足误差限。

MATLAB求解方程解析解和数值解

MATLAB求解方程解析解和数值解

辽宁工程技术大学上机实验报告用MATLAB求解质点振动方程振动是日常生活和工程技术中常见的一种运动形式。

利用常系数线性微分方程的理论来讨论有关自由振动和强迫振动的相关问题。

利用MA TLAB数学软件大致可分四类情况:(1)无阻尼自由振动情况;(2)有阻尼自由振动;(3)无阻尼强迫振动;(4)有阻尼强迫振动求其数值解和解析解;MATLAB软件求解微分方程解析解的命令“dsolve()”求通解的命令格式:(’微分方程’,’自变量’)注:微分方程在输入时,一阶导数y’应输入Dy,y’’应输入D2y等,D应大写。

1,无阻尼自由振动情况:常见的数学摆的无阻尼微小振动方程代码如下:>> t=0:pi/50:2*pi;>> y=2*sin(3*t+2);>>plot(t,y,'b')2,有阻尼自由振动由无阻尼振动的通解可以看出,无阻尼振动是按照正弦规律运动的,摆动似乎可以无限期的进行下去,但事实上,空气从在阻力,在运动时,我们必须把空气阻力考虑在内,所以我们得到有阻尼摆动方程为:记u/m=2n,g/l=w^2,这里n,w是正常数,所以:y=dsolve('D2y+2*n*Dy+w^2*y=0','t'); (4.43)解得:y = C3*exp(-t*(n + ((n + w)*(n - w))^(1/2))) + C2*exp(-t*(n - ((n + w)*(n - w))^(1/2)))(1)小阻尼情形:n<w时,方程(4.43)的通解为:y=exp(-n*t)*(c1*cos(w1*t)+c2*sin(w1*t))和前面无阻尼的情形一样,可以把上式的通解改写为一下形式:y=A*exp(-n*t)*sin(w1*t+Q), (4.45)这里的A,Q为任意常数。

用matlab 操作得到:t=0:0.1:10;y=3*exp(-0.1*t).*sin(5*t+4);plot(t,y,'k-')如图:由(4.45)可见,摆动的运动不是周期的,振动的幅度随着时间的增加而不断减小。

MATLAB方程组求解实验报告

MATLAB方程组求解实验报告
通过实验,学会了用MATLAB语言编写程序求解线性、非线性方程组的解。了解了线性、非线性方程组用编程求解方法,运用MATLAB语言编写应用程序,完成对线性、非线性方程组的求解更加方便快捷。
八、教师评语
签名:
日期: 年 月 日
成绩
14 -2 2
0 33 8
-2 7 38
3)
a=[2 1 1;1 3 1;1 1 4];
b=[12 -3 1;-1 30 7;-3 6 34];
>> a-b
ans =
-10 4 0
2 -27 -6
4 -5 -30
4)
a=[2 1 1;1 3 1;1 1 4];
b=[12 -3 1;-1 30 7;-3 6 34];
方程A*X=b变形成QRX=b则X=R\(Q\b)
命令[Q,R]=qr(A)
2求线性齐次方程组的通解
在Matlab中,函数null用来求解零空间,即满足A·X=0的解空间,实际上是求出解空间的一组基(基础解系)。
格式:z=null %z的列向量为方程组的正交规范基,满足
% z的列向量是方程AX=0的有理基
1)
a=[2 1 1;1 3 1;1 1 4];
b=[12 -3 1;-1 30 7;-3 6 34];
>> a*b
ans =
20 30 43
6 93 56
-1 51 144
2)
a=[2 1 1;1 3 1;1 1 4];
b=[12 -3 1;-1 30 7;-3 6 34];
>> a+b
ans =
-2.3529 12.5294 0.7059
-2.2353 0.3529 9.4706

数值分析matlab实验报告

数值分析matlab实验报告

数值分析matlab实验报告《数值分析MATLAB实验报告》摘要:本实验报告基于MATLAB软件进行了数值分析实验,通过对不同数学问题的数值计算和分析,验证了数值分析方法的有效性和准确性。

实验结果表明,MATLAB在数值分析领域具有较高的应用价值和实用性。

一、引言数值分析是一门研究利用计算机进行数值计算和分析的学科,其应用范围涵盖了数学、物理、工程等多个领域。

MATLAB是一种常用的数值计算软件,具有强大的数值分析功能,能够进行高效、准确的数值计算和分析,因此在科学研究和工程实践中得到了广泛的应用。

二、实验目的本实验旨在通过MATLAB软件对数值分析方法进行实验验证,探究其在不同数学问题上的应用效果和准确性,为数值分析方法的实际应用提供参考和指导。

三、实验内容1. 利用MATLAB进行方程求解实验在该实验中,利用MATLAB对给定的方程进行求解,比较数值解和解析解的差异,验证数值解的准确性和可靠性。

2. 利用MATLAB进行数值积分实验通过MATLAB对给定函数进行数值积分,比较数值积分结果和解析积分结果,验证数值积分的精度和稳定性。

3. 利用MATLAB进行常微分方程数值解实验通过MATLAB对给定的常微分方程进行数值解,比较数值解和解析解的差异,验证数值解的准确性和可靠性。

四、实验结果与分析通过对以上实验内容的实际操作和分析,得出以下结论:1. 在方程求解实验中,MATLAB给出的数值解与解析解基本吻合,验证了MATLAB在方程求解方面的高准确性和可靠性。

2. 在数值积分实验中,MATLAB给出的数值积分结果与解析积分结果基本吻合,验证了MATLAB在数值积分方面的高精度和稳定性。

3. 在常微分方程数值解实验中,MATLAB给出的数值解与解析解基本吻合,验证了MATLAB在常微分方程数值解方面的高准确性和可靠性。

五、结论与展望本实验通过MATLAB软件对数值分析方法进行了实验验证,得出了数值分析方法在不同数学问题上的高准确性和可靠性。

数值分析实验报告matlab

数值分析实验报告matlab

数值分析实验报告matlab数值分析实验报告引言:数值分析是一门研究利用计算机数值方法解决数学问题的学科,它在科学计算、工程设计、金融分析等领域具有重要的应用价值。

本实验报告旨在通过使用MATLAB软件,探索数值分析的基本原理和方法,并通过实际案例加深对数值分析的理解。

一、误差分析在数值计算中,误差是无法避免的。

误差分析是数值分析中的重要一环,它帮助我们了解数值计算的准确性和稳定性。

在实验中,我们通过计算机模拟了一个简单的数学问题,并分别计算了绝对误差和相对误差。

通过比较不同算法的误差大小,我们可以选择最适合的算法来解决实际问题。

二、插值与拟合插值和拟合是数值分析中常用的方法,它们可以通过已知的数据点来推导出未知数据点的近似值。

在本实验中,我们通过MATLAB的插值函数和拟合函数,分别进行了插值和拟合的实验。

通过比较不同插值和拟合方法的结果,我们可以选择最适合的方法来处理实际问题。

三、数值积分数值积分是数值分析中的重要内容,它可以用来计算曲线下的面积或函数的积分值。

在实验中,我们通过MATLAB的数值积分函数,对一些简单的函数进行了积分计算。

通过比较数值积分和解析积分的结果,我们可以评估数值积分的准确性和稳定性,并选择最适合的积分方法来解决实际问题。

四、常微分方程的数值解法常微分方程是数值分析中的重要内容,它可以用来描述许多自然现象和工程问题。

在实验中,我们通过MATLAB的常微分方程求解函数,对一些简单的微分方程进行了数值解法的计算。

通过比较数值解和解析解的结果,我们可以评估数值解法的准确性和稳定性,并选择最适合的数值解法来解决实际问题。

五、线性方程组的数值解法线性方程组是数值分析中的经典问题,它在科学计算和工程设计中广泛应用。

在实验中,我们通过MATLAB的线性方程组求解函数,对一些简单的线性方程组进行了数值解法的计算。

通过比较数值解和解析解的结果,我们可以评估数值解法的准确性和稳定性,并选择最适合的数值解法来解决实际问题。

  1. 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
  2. 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
  3. 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。

佛山科学技术学院
实 验 报 告
课程名称 数值分析
实验项目 数值积分
专业班级 机械工程 姓 名 余红杰 学 号 2111505010 指导教师 陈剑 成 绩 日 期 月 日
一、实验目的
1、 掌握程序的录入和matlab 的使用和操作;
2、 了解影响线性方程组解的精度的因素——方法与问题的性态。

3、 学会Matlab 提供的“\”的求解线性方程组。

二、实验要求
1、按照题目要求完成实验内容;
2、写出相应的Matlab 程序;
3、给出实验结果(可以用表格展示实验结果);
4、分析和讨论实验结果并提出可能的优化实验。

5、写出实验报告。

三、实验步骤
1、用LU 分解及列主元高斯消去法解线性方程组
a)⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎪⎪⎭⎫ ⎝⎛----15900001.582012151526099999.2310
7104321x x x x , 输出b Ax =中系数LU A =分解的矩阵L 和U ,解向量x 和)det(A ;用列主元法的行交换次序解向量x 和求)det(A ;比较两种方法所得结果。

2、用列主高斯消元法解线性方程组b Ax =。

(1)、⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4987.023.116.427
.199.103.601.3321x x x
(2)、⎪⎪⎪⎭
⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎪⎪⎪⎭⎫ ⎝⎛--11134.981.4990.023.116.427
.199.103.600.3321x x x 分别输出)det(,,A b A ,解向量x ,(1)中A 的条件数。

分析比较(1)、(2)的计算结果
3、线性方程组b Ax =的A 和b 分别为
⎪⎪⎪⎪⎪⎭⎫ ⎝⎛=1095791068565778710A ,⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=31332332b 则解T x ),1,1,1,1(=. 用MATLAB 内部函数求)det(A 和A 的所有特征值和2)(A cond . 若令
⎪⎪⎪⎪⎪⎭
⎫ ⎝⎛=+98.99599.6989.998.585604.508.72.71.8710A A δ, 求解b x x A A =++))((δδ,输出向量x δ和2x δ,从理论结果和实际计算两方面分析线性
方程组b Ax =解的相对误差22/x x
δ以及A 的相对误差
/A A δ的关系。

四、实验结果
1:
%run311.m
clc,clear;
A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];
b = [8;5.90001;5;1];
%L U 分解
format short
%小数点后四位,不然会受到后面的影响
[L U] = lu(A)
%解方程组,输出A ,det(A)
y = L\b;
format long
%小数点后15位显示
x = U\y
%V ol_xiao.m列主元消去,未用增广矩阵而把系数矩阵A和矩阵B分开对应来的function x=V ol_xiao(A,b)
%列主元消去
%x为解
%A为系数矩阵
n =length(A);
%A默认为方阵,求其大小,且默认和b长度一样
C =zeros(1,n);
%用于系数矩阵交换存放数据
b0 =0;
%用于b矩阵交换存放数据
x =zeros(n,1);
%建一个矩阵存放解
%下面对矩阵进行n-1次列主元交换得到上三角矩阵
for i=1:(n-1)
for j=(i+1):n
if abs(A(i,i))<abs(A(j,i))
C =A(i,1:n);
A(i,1:n)=A(j,1:n);
A(j,1:n)=C;
b0 = b(i);
b(i) =b(j);
b(j) =b0;
else continue;
end
end;
for j=(i+1):n
b(j) = b(j) + b(i)*(-A(j,i))/A(i,i);
A(j,1:n)=A(j,1:n)+A(i,1:n).*(-A(j,i)/A(i,i));
end
end
%以下为判断是否满秩然后从下往上求解:
if A(n,n)~=0
x(n) = b(n)/A(n,n);
for i=(n-1):-1:1
sum = 0;
for j=(i+1):n
sum = sum + A(i,j)*x(j);
end
x(i) = (b(i)-sum)/A(i,i);
end
else x = ['err'];
end
%run312.m
clc,clear;
A = [10 -7 0 1;-3 2.099999 6 2;5 -1 5 -1;2 1 0 2];
b = [8;5.90001;5;1];
format long
x=V ol_xiao(A,b)
(其实可以看出,两种结果是一致的)
2.
%run321.m
clc;clear;
A = [3.01 6.03 1.99;1.27 4.16 -1.23;0.987 -4.81 9.34];
b = [1;1;1];
format bank
%格式转换为普通的,不然会用科学计数法表示
x = Vol_xiao(A,b)
T1 = cond(A,1)
T2 = cond(A,2)
T3 = cond(A,inf)
%可以看出条件数很大
%run322.m
clc;clear;
A = [3.01 6.03 1.99;1.27 4.16 -1.23;0.990 -4.81 9.34]
b = [1;1;1]
det1 = det(A)
format bank
x = Vol_xiao(A,b)
%从条件数就可以猜出其实微小变化就会对结果造成很大差异
%结果也验证了条件数大的方程的不稳定性很差。

3.
%run331.m
clc;clear;
A = [10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
b = [32;23;33;31];
format short
%显示对应行列式值和特征值
det1 = det(A)
eig1 = eig(A)
format bank
%显示条件数
tjs2 = cond(A,2)
%run332.m
clc;clear;
A = [10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10];
b = [32;23;33;31];
X = [1;1;1;1];
AA= [10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 5 9 9.98];
oA= AA - A;
XX= AA \ b;
oX= XX – X
%输出oX 的2范数
normx2 = norm(oX)
%输出X 的相对误差
xrerro = normx2/norm(X)
%输出A 的相对误差
Arerro = norm(oA)/norm(A)
五、讨论分析
1.对于第一题,有一个顺序主子式是极小的,接近等于0,所以无法直接得到标准的上下三角,所以调用程序内置的lu 函数进行分解,所求到的值与列主元方法的值一致,查找后发现matlab 的lu 分解就是根据列主元方法求出的;
2.
条件数对一个方程的稳定性有很大影响,题中的系数矩阵条件数超过五万,微
小的变化就引起了结果的很大变化,如果是数据分析的话相当于蝴蝶效应般严重;
3.第三题也是对于一个系数矩阵的范数,很小的变化引起结果的大变化,该系数矩阵的相关数也是三千多,也是比较大,不稳定的。

六、改进实验建议
相关题目太少,而且很少有扩展性题目,对于各种指令的使用并不充分,其实在这个求解方程组的过程中,还是有比较多的相关方法可以进行MATLAB实现的。

相关文档
最新文档