当前位置：文档之家 › 大连理工大学矩阵与数值分析MATLAB上机实验

大连理工大学矩阵与数值分析MATLAB上机实验

  • 格式：pdf
  • 大小：912.87 KB
  • 文档页数：18

二、解线性方程组 1．分别 Jacobi 迭代法和 Gauss-Seidel 迭代法求解线性方程组
3 1 0 0 x1 1 1 3 1 0 x2 0 , 0 1 2 1 x3 0 0 0 1 3 x4 0
Gauss 列主元消去法程序:
clc; clear; format long a=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; %系数矩阵 b=[12;6;23;16]; [n,m]=size(a); nb=length(b); det=1; for k=1:n-1 amax=0; for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k)); r=i; end end if amax<1e-10 return; end if r>k for j=k:n z=a(k,j); a(k,j)=a(r,j); a(r,j)=z; end z=b(k);
从小到大求和程序计算结果：
N 100 10000 1000000 从小到大求和程序得 到的 SN 0.497512437810945 0.499975001249937 0.499999750000134 真实值������������ = ������ 2������ + 1 0.497512437810945 0.499975001249937 0.499999750000125 计算值有效位数 15 15 13
8
2
1 dx x
复化梯形公式程序
clc; clear; format long syms t m=int(1/t,2,8); %真实值 a=2; b=8; n=300; h=(b-a)/n; sum=0; f=inline('1/x'); for i=1:n-1 sum=sum+f(a+i.*h); end T=h/2*(f(a)+2*sum+f(b))
b(k)=b(r); b(r)=z; det=-det; end for i=k+1:n m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); end det=det*a(n,n); for k=n:-1:1 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end x=x'
f ' ( xk ) f ( xk ) f ( xk 1 ) ，即用割线代替切线，这引入了误差。 xk xk 1
2、生成一个列和为 1 的 100 阶随机矩阵，编写幂法程序求其最大特征值。
幂法程序
clc; clear; format long A=rand(100,100); for i=1:100 A(:,i)=A(:,i)/sum(A(:,i)); end N=100; %最大迭代次数 ep=1e-6; %精度要求 n=length(A); u=ones(n,1); k=0; m1=0; while k<=N v=A*u; m=max(abs(v)); u=v/m; if abs(m-m1)<ep break; end m1=m; k=k+1; %生成100阶随机矩阵 %使随机矩阵；列和为1
一、设 S N
j 1
N
1 4 j2 1
，分别编制从小到大和从大到小的顺序程序计算
S100 , S10000 , S1000000 ，
并指出有效位数。 从小到大求和程序：
clc; clear; format long S=0; N=100; %此处N分别取100,10000,1000000进行计算 for j=1:N S=S+1/(4*j^2-1); end S
复化梯形公式程序计算结果
分析
本程序将区间[2,8]进行了 300 等分，在命令行输入 abs(T-m)，可得计算值与精确 值的误差为 7.812416996655358e-06，误差未超过 10-5，满足题目要求。
三点 Gauss 型公式程序
clc; clear; fun=inline('1./x'); a = 2; b = 8; tol=1e-5; n=2; syms x %精度 %n次代数精度，n=(3-1)*2+1 %被积函数
的根，计算停止的条件为： xk 1 xk 106
Newton 迭代法程序
clc; clear; format long x0=1; %设初始值
f=inline('exp(x)+2^(-x)+2*cos(x)-6'); diff=inline('exp(x)-2^(-x)*log(2)-2*sin(x)'); x=x0-feval(f,x0)/feval(diff,x0); while abs(x-x0)>=1e-6 x0=x; x=x0-(exp(x0)+2^(-x0)+2*cos(x0)-6)/(exp(x0)-2^(-x0)*log(2)-2*sin(x0)) ; end x
矩阵与数值分析之 MATLAB 试验
学生姓名：陈贵鹏 学 号：21606075
专业班级：结构工程建研 1602 班 所在单位：智能结构研究所
2017 年 6 月 4 日ຫໍສະໝຸດ 2016 级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
教学班号：建研 1602 班 任课教师：董波 学生姓名：陈贵鹏 学号：21606075 院系：建设工程学部结构工程系
迭代法计算停止的条件为： max x (jk 1) x (jk ) 10 6 ．
1 j 3
Jacobi 迭代程序:
clc; clear; format long A=[-3,1,0,0;1,-3,1,0;0,1,-2,1;0,0,1,-3]; b=[-1;0;0;0]; D=diag(diag(A)); %求对角矩阵D L=-tril(A,-1); U=-triu(A,1); B=D\(L+U); f=D\b; x0=[0;0;0;0]; x=B*x0+f; x0=x; %取零向量作为初始向量 %求下三角矩阵L %求上三角矩阵U
QR 方法程序运算结果：
分析
从计算结果来看，如果不计舍入误差，Gauss 列主元消去法与 QR 方法得到的结 果几乎相同。但是从编程的复杂程度来看，个人感觉 Gauss 列主元消去法难度更 大，QR 方法的语句相对简洁。
三、非线性方程的迭代解法
1、用 Newton 迭代法、割线法求方程
f x e x 2 x 2cos x 6 0
Newton 迭代法程序运算结果
割线法程序
clc; clear; format long xm=1; %设初始值 xn=0.5; %设初始值 f=inline('exp(x)+2^(-x)+2*cos(x)-6'); while abs(xn-xm)>=1e-6 x=xn-(feval(f,xn)*(xn-xm))/(feval(f,xn)-feval(f,xm)); xm=xn; xn=x; end x
while norm(x-x0,'inf')>=1e-6
x=B*x0+f; end x
Jacobi 迭代程序运算结果：
Gauss-Seidel 迭代程序:
clc; clear; format long A=[-3,1,0,0;1,-3,1,0;0,1,-2,1;0,0,1,-3]; b=[-1;0;0;0]; D=diag(diag(A)); %求A的对角矩阵D L=-tril(A,-1); U=-triu(A,1); G=(D-L)\U; f=(D-L)\b; x0=[0;0;0;0]; x=G*x0+f; x0=x; x=G*x0+f; end x %取零向量作为初始向量 %求下三角矩阵L %求上三角矩阵U
割线法程序运算结果
分析
在 MATLAB 命令窗口输入 x=solve('exp(x)+2^(-x)+2*cos(x)-6=0')，可以得到方程 的根为 x= 1.8293836019338488171362129468142。对比 Newton 法和割线法求得 的结果，可以发现 Newton 法求得的结果更精确。这是因为割线法采用了导数的 近似式
2. 用 Gauss 列主元消去法、QR 方法求解如下方程组：
2 8 5 9 4 3 1 x1 12 2 0 0 x2 6 , 0 4 0 x3 23 0 0 5 x4 16
分析：
由 SN
1 1 1 N = ( ) 可以获得真实值。当 N 很大时， 2 j 1 2N 1 j 1 4 j 1 j 1 2 2 j 1
2
N
1
N
可以发现从小到大计算与从大到小的计算结果出现了差异。 这是因为从大到小求 和时，变量在分母上，整个算式是由小到大进行求和的，这样求和避免了“大数 吃小数”现象，因而计算结果更优，具有更高的有效位数。
Gauss 列主元消去法程序运算结果:
QR 方法程序：
clc; clear; format long A=[2,4,3,1;8,2,0,0;5,0,4,0;9,0,0,5]; %系数矩阵 B=[12;6;23;16]; n=length(B); Q=eye(4); H=A; for i=1:n-1 b=H(:,1);
c=norm(b,2); b(1)=b(1)-c; b1=transpose(b); Q1=eye(4-i+1)-2/(b1*b)*b*b1; H1=Q1*H; Q1=blkdiag(eye(i-1),Q1); Q=Q*transpose(Q1); if i<n-1 H1(1,:)=[]; H1(:,1)=[]; H=H1; end end R=transpose(Q)*A; x=R\(Q\B)

合集下载

大连理工矩阵上机作业

大连理工矩阵上机作业

第一题Lagrange插值函数function y=lagrange(x0,y0,x);n=length(x0);m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif j~=kp=p*(z-x0(j))/(x0(k)-x0(j));endends=p*y0(k)+s;endy(i)=s;endx0=[1:10];y0=[67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76 .777];lagrange(x0,y0,17)ans= -1.9516e+12x=[1:0.1:10];x=x';plot(x0,y0,'r');hold onplot(x,y,'k');legend('原函数','拟合函数')拟合图像如下拟合函数出现了龙格现象，运用多项式进行插值拟合时，效果并不好，高次多项式会因为误差的不断累积，导致龙格现象的发生。

第二题function fun =nihe(n)m=[67.052*10^6,68.008*10^6,69.803*10^6,72.024*10^6,73.400*10^6,72.063 *10^6,74.669*10^6,74.487*10^6,74.065*10^6,76.777*10^6];w=[1,2,3,4,5,6,7,8,9,10];d1=0;d2=0;d3=0;y1=polyfit(m,w,1);y2=polyfit(m,w,2);y3=polyfit(m,w,3);y2=poly2sym(s2);y3=poly2sym(s3);y4=poly2sym(s4);f1=subs(y1,17);f2=subs(y2,17);f3=subs(y3,17);for p=1:10;d1=d1+(subs(y1,w(p))-m(p))^2;d2=d2+(subs(y2,w(p))-m(p))^2;d3=d3+(subs(y3,w(p))-m(p))^2;endd1=sqrt(d1);d2=sqrt(d2);d3=sqrt(d3);fun=[f1 f2 f3;d2 d3 d4];return;结果拟合形式预计2010年产量均方误差直线83382272.733087388.906抛物线74410590.912601531.997三次曲线99041977.162396531.259三次函数的均方误差最小，拟合的最好。

MATLAB实验上机所用word(数值计算篇)

MATLAB实验上机所用word(数值计算篇)

第三章数值计算3.1LU分解和恰定方程组的解Matlab矩阵的分解形式主要有：三角分解、正交化、特征值分解。

3.1.1LU三角分解、行列式和逆1、LU分解是利用Gaussian列主元消去法进行的。

为了保证主元消去策略的实施，一般来说，必须对被分解矩阵实施行置换因此有：LU=PA 式中L为主对角元为1的下三角矩阵，U是上三角矩阵，P是由0或1组成的行置换矩阵：［L UP］=lu(A)2、A-1=U-1L-1P Matlab根据以上算法编制了相应的指令：det(a) 求矩阵a的行列式inv(a) 求矩阵a的逆矩阵3.1.2恰定方程组的解在求解方程式，尽量不要使用inv(a)*b指令，而应采用a\b，后者计算速度比前者快，精确度高。

【例3.1.2-1】“求逆”法和“左除”法解恰定方程的性能对比（1）randn('state',0);A=gallery('randsvd',100,2e13,2); ％产生条件数为2e13的100阶随机矩阵x=ones(100,1);b=A*x;（2）ticxi=inv(A)*b;ti=toc（3）tic;xd=A\b;td=toc;3.1.3矩阵特征值和矩阵函数矩阵A与向量x相乘,即表示矩阵对向量的变换(transformation)，一般说来，向量在变换的作用下将发生旋转（rotation），反射(reflection)和放大缩小。

但对于任何一个矩阵来说，中存在那么一些特殊的向量，再对其变化的作用下，向量的方向不变，而仅长短发生变化。

这种向量就是所谓的特征向量（eigenvector），它满足方程 Ax=λx3.1.3.1特征值和特征向量的求取d=eig(a) 仅仅计算a的特征值[v d]=eig(a) 计算矩阵a的特征向量阵v和特征对角阵d，使av=vd成立。

【例3.1.3.1-1】简单实阵的特征值问题。

A=[1,-3;2,2/3];[V,D]=eig(A)V =0.7746 0.77460.0430 - 0.6310i 0.0430 + 0.6310iD =0.8333 + 2.4438i 00 0.8333 - 2.4438i3.2数据分析3.2.1 基本统计函数指令note: 1、median(x) 当x为向量时，先把x元素有小到大排列，然后在新排成数组中；n为奇取（n+1）/2个元素；n为偶数时取n/2与(n/2+1)元素的平均值，作为总位数。

大连理工大学工科数学分析上机作业

大连理工大学工科数学分析上机作业

工科数学分析上机作业说明：以下两道题均是使用Matlab 语言，且在Matlab 7.0中运行通过。

1.（两个重要极限）计算下列函数的函数值并画出图形，观察两个重要极限值。

(1)y=f(x)=; (2)y=f(x)=.解：（1）求解过程如下：>> syms x>> y=limit(sin(x)/x)y =1>> ezplot(sin(x)/x,[-10*pi,10*pi])>> ezplot(sin(x)/x,[-1*pi,1*pi])其图形如下：（2）求解过程如下：>> syms x>> y=(1+x)^(1/x)y =(1+x)^(1/x)>> y=limit((1+x)^(1/x))y =exp(1)>> ezplot((1+x)^(1/x),[-1000,1000]) >> ezplot((1+x)^(1/x),[-10,10]) >> ezplot((1+x)^(1/x),[-1,1])其图像如下：分析如下：（1）当x 取值为[-30,30]时，由该题的第一个图像可以看到，函数值在不断震荡，一会为正数，一会为负数。

而当x取值为[-3，3]时，函数值始终大于0。

当x趋近于0时，由该题的第二个图像可以得到函数值为1。

另外，该结论也可以由夹逼法则证明，结果不变，当x趋近于0时，函数值仍为1。

（2）由该题的三个图像可以知道，该函数在定义域内为单调递减函数。

且由该题的第一和二个图像知道，当x在[0,10]区间内，函数递减趋势非常迅速。

由该题的第三个图像知道，当x趋于0 时，函数值为自然对数的底数 e ,即约为2.71828.3.计算f(x)=,的函数值{f（0.1k）;k=1,2,…,30}.计算结果取7位有效数字。

解：计算过程为：>> f1= @(t) exp(-(t).^2/2)f1 =@(t) exp(-(t).^2/2)>> for i=1:30s(i,1)=1/2+1/sqrt(2*pi)*quad(f1,0,i*0.1);end>> fprintf('%9.7g\n',s);0.53982780.57925970.61791140.65542170.69146250.72574690.75803630.78814460.81593990.84134470.86433390.88493030.90319950.91924340.93319290.94520070.95543450.96406970.97128340.97724990.98213560.98609660.98927590.99180250.99379030.99533880.99653290.99744480.99813420.9986501>>分析：本题使用了Matlab 强大的数学计算（微积分计算）功能。

MATLAB矩阵实验报告

MATLAB矩阵实验报告

MATLAB程序设计实验班级：电信1104班姓名：龙刚学号：1404110427实验内容：了解MA TLAB基本使用方法和矩阵的操作一．实验目的1.了解MA TLAB的基本使用方法。

2.掌握MA TLAB数据对象的特点和运算规则。

3.掌握MA TLAB中建立矩阵的方法和矩阵的处理方法。

二．实验内容1.浏览MATLAB的start菜单，了解所安装的模块和功能。

2.建立自己的工作目录，使用MA TLAB将其设置为当前工作目录。

使用path命令和工作区浏览两种方法。

3.使用Help帮助功能，查询inv、plot、max、round等函数的用法和功能。

使用help命令和help菜单。

4.建立一组变量，如x=0:pi/10:2*pi，y=sin(x)，在命令窗口显示这些变量；在变量窗口打开这些变量，观察其值并使用绘图菜单绘制y。

5.分多行输入一个MA TLAB命令。

6.求表达式的值)610.3424510w-=+⨯()22tanb ca eabcxb c aππ++-+=++，a=3.5，b=5，c=-9.8(20.5ln tz e t=，21350.65it-⎡⎤=⎢⎥-⎣⎦7.已知1540783617A--⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦，831253320B-⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦求A+6B，A2-B+IA*B，A.*B，B*AA/B，B/A[A,B]，[A([1,3], :); B^2]8.已知23100.7780414565532503269.5454 3.14A -⎡⎤⎢⎥-⎢⎥=⎢⎥⎢⎥-⎣⎦ 输出A 在[10,25]范围内的全部元素取出A 的前三行构成矩阵B ，前两列构成矩阵C ，右下角3x2子矩阵构成矩阵D ，B 与C 的乘积构成矩阵E分别求表达式E<D ，E&D ，E|D ，(~E) | (~D)9.已知2961820512885A -⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦求A 的特征值和特征向量，分析其数学意义。