当前位置：文档之家 › 大连理工大学矩阵与数值分析上机作业代码

大连理工大学矩阵与数值分析上机作业代码

  • 格式：pdf
  • 大小：372.27 KB
  • 文档页数：37

n= 9 （3）取 x
(0)
= 0 ， b =（1,0,1,0,……，1,0）T,迭代误差取默认值
10
−5
，迭代次数上限取默认
值 50，使用 Jacobi 法进行迭代。 >> test2 >> b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]' >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.3238 -0.0985 0.3116 -0.0881 0.3110
…… 1.0000 1.0000 1.0000 1.0000
1.2 计算结果分析
对于高斯消元法， 当 n=10 时， 得到的计算结果为 x= （1.0000， 1.0000， 1.0000， 1.0000， ……， T 1.0000） ， 恰好为精确值； 但当 n=84 时， 得到的计算结果却为 x=(0.0000,0.0000,0.0000,……， -0.0001,0.0002， ……，2.3487，-4.0263，5.3684）T,完全偏离了精确值。 对于列主元消去法，当 n=10 时，得到的计算结果为 x= （ 1.0000 ， 1.0000 ， 1.0000 ， T 1.0000， ……， 1.0000） ， 恰好为精确值； 当 n=10 时， 得到的计算结果为 x= （1.0000， 1.0000， 1.0000，1.0000， ……，1.0000）T，恰好为精确值。 由于使用高斯消元法时可能会出现小主元作除数的情况，虽然在求解系数矩阵为 10 阶 的方程时，由于步数较少，尚能求得精确解，但在求解求解系数矩阵为 84 阶的巨型方程时， 由于求解步数太多，使得小主元作除数的情况累加，最终导致“大数除小数”的情况发生， 因而所得解严重偏离精确解。 而使用列主元消去法时， 由于有效地避免了这一不利情况， 因而无论是在求解系数矩阵 为 10 阶还是 84 阶的方程时，都能得到精确解。
-0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0882 0.3105 -0.0957 0.3313 -0.0472
n= 16 （4）取 x
(0)
= 0 ， b =（1,0,1,0,……，1,0）T,迭代误差取默认值
习题 2
设有方程组 Ax = b ，其中 A ∈ R
20× 20
，
−0.5 −0.25 ⎛ 3 ⎜ 3 −0.5 ⎜ −0.5 ⎜ −0.25 −0.5 ⋱ ⎜ A=⎜ ⋱ ⋱ ⎜ ⋱ ⎜ ⎜ ⎜ ⎝
(a) 选 取 不 同 的 初 始 向 量 x
(0)
⋱ ⋱ ⋱ ⋱ ⋱
⎞ ⎟ ⎟ ⎟ ⋱ ⎟ ⋱ ⎟ ⋱ −0.5 −0.25 ⎟ ⎟ −0.5 3 −0.5 ⎟ −0.25 −0.5 3 ⎟ ⎠
-0.0209 0.0419 -0.0836 0.1665 -0.3303 0.6501 -1.2582 2.3487 -4.0263 5.3684 （2）列主元消去法 n=10 时， >> [A,b]=CreateMatrix(10) >> x=MainElementMethod(A,b) x= 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 n=84 时， >> [A,b]=CreateMatrix(84) >> x=MainElementMethod(A,b) x= 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
U=p{i}*U; for j=i+1:n l{i}(j,i)=-U(j,i)/Max; end U=l{i}*U; end %计算L L=eye(n); for i=1:n-2 for j=i+1:n-1 l{i}=p{j}*l{i}/p{j}; end L=l{i}*L; end L=inv(L*l{n-1}); %计算P P=eye(n); for i=1:n-1 P=p{i}*P; end %计算x x=U\(L\P*b);
1.3.2 GaussMethod.m function x=GaussMethod(A,b) %高斯消元法返回方程Ax=b的解
%x是方程的解 %A为方程的系数矩阵 %b为方程的常数矩阵 n=length(b); %定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵 l=cell(1,n-1); l(:)={eye(n)}; U=A; for i=1:n-1 for j=i+1:n l{i}(j,i)=-U(j,i)/U(i,i); end U=l{i}*U; end %计算L和U L=inv(l{1}); for i=2:n-1 L=L/l{i}; end x=U\(L\b); .m 1.3.3 MainElementMethod MainElementMethod.m function x=MainElementMethod(A,b) %列主元消去法返回方程Ax=b的解 %x是方程的解 %A为方程的系数矩阵 %b为方程的常数矩阵 n=length(b); %定义单位下三角矩阵元胞数组，并将全部数组元素赋值为n阶单位阵 l=cell(1,n-1); l(:)={eye(n)}; p=cell(1,n-1); p(:)={eye(n)}; U=A; for i=1:n-1 %找出列主元所在行，并交换首行，得到Pi [Max,row]=max(U(i:n,i)); row=row+i-1; mid=p{i}(row,:); p{i}(row,:)=p{i}(i,:); p{i}(i,:)=mid;
1.3 M 文件
.m 1.3.1 CreateMatrix CreateMatrix.m function [A,b]=CreateMatrix(n) %用于存放习题1的题目信息，并建立构造题目中矩阵的函数 %对矩阵A赋值 A1=6*ones(1,n); A2=ones(1,n-1); A3=8*ones(1,n-1); A=diag(A1)+diag(A2,1)+diag(A3,-1); %对向量b赋值 b=15*ones(n,1); b(1)=7; b(n)=14;
0.3111 -0.0882 0.3104 -0.0957 0.3313 -0.0472 n= 9 （5）取 x
(0)
=（1,1， ……，1,1）T， b =（1,1， ……，1,1）T,迭代误差取默认值
10
−5
，迭代次
数上限取默认值 50，使用 Jacobi 法进行迭代。 >> test2 >> b=zeros(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2228 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2228 0.2223 0.2159 0.2327 0.2766 n= 19 （6）取 x
将 A 的主对角线元素成倍增长若干次，非主对角元素不变，每次用 Jacobi 法计算，要求迭 代误差满足 x ( k +1) − x ( k )
∞
< 10−6 ，比较收敛速度，分析现象并得出你的结论。
2.1 程序
(0) −5
（1）取 x
= 0 ， b =（1,1， ……，1,1）T,迭代误差取默认值
10
−5
，迭代次数上限取默认
值 50，使用 Gauss-Seidel 法进行迭代。 >> test2 >> b=[1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0]' >> x0=zeros(20,1) >> [x,n]=Gauss_SeidelMethod(A,b,x0) x= 0.3238 -0.0985 0.3115 -0.0881 0.3110 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889 0.3111 -0.0889
n= 17 （2）取 x
(0)
= 0 ， b =（1,1， ……，1,1）T,迭代误差取默认值
10
−5
，迭代次数上限取默认值
50，使用 Gauss-Seidel 法进行迭代。 >> test2
>> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=Gauss_SeidelMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2228 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2228 0.2223 0.2159 0.2327 0.2766
10
，迭代次数上限取默认值
50，使用 Jacobi 法进行迭代。 >> test2 >> b=ones(20,1) >> x0=zeros(20,1) >> [x,n]=JacobiMethod(A,b,x0) x= 0.2766 0.2327 0.2159 0.2223 0.2227 0.2221 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2222 0.2221 0.2227 0.2223 0.2159 0.2327 0.2766

合集下载

大连理工线代上机实验

大连理工线代上机实验


基本操作
四则运算、转置、求逆、求秩、求行列式、组合、 化为行最简形、求特征值
常见任务
① 矩阵的赋值和其加、减、乘、除（求逆）命令； ② 矩阵化为最简行阶梯型的计算命令；[U0,ip]=rref(A) ③ 多元线性方程组MATLAB求解的几种方法；x=inv(A)*b, U=rref(A) ④ 行列式的几种计算机求解方法； D=det(A),[L,U]=lu(A);D=prod(diag(L)) ⑤ n个m维向量组的相关性及其秩的计算方法和命令; r=rank(A),U=rref(A) ⑥ 求欠定线性方程组的基础解系及超定方程解的MATLAB 命令；xb=null(A) ⑦ 矩阵的特征方程、特征根和特征向量的计算命令； f=poly(A);[P,D]=eig(A) ⑧ 化二次型为标准型的MATLAB命令；yTDy=xTAx; 其中 y=P-1x,
• • • •
-66.5556 25.6667 -18.7778 26.5556
• >>
例三、求秩
• • • • • • • • • • • • >> A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; >> r=rank(A); >> r= % = 计算机不显示r的值 ??? r= | Error: Expression or statement is incomplete or incorrect. >> rank(A) ans = 4 >> r r= 4 %不打；则计算机将显示rank(A)的值
例1 用直接解法求解下列线性方程组. 命令如下： A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4]; b=[13,-9,6,0]'; x=A\b
Removed_大连理工大学工科数学分析上机作业

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

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

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

(1)y=f(x)=; (2)y=f(x)=.sin x x (1+x)1x 解：（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)=,12+1√2π∫x 0e ‒t 2/2dt 1≪x ≪3的函数值{f （0.1k ）;k=1,2,…,30}.计算结果取7位有效数字。

矩阵与数值分析 误差分析

矩阵与数值分析 误差分析


x
称为近似值a 的相对误差。 相对误差也可正可负。 实际计算中，如果真值 x 未知时， 通常取
作为a的相对误差， 条件是
xa xa x a
xa 较小。 x
相对误差的绝对值上界叫做相对误差界（限）， 记为：
xa a
ea a
相对误差界（限）
例
有两个量 x=3.000， a=3.100， 则其绝对误差： x a 0.1 绝对误差 其相对误差为： x a 0.1 0.333 101 , x 3.00 又有两个量 x 300.0 ， a 310.0, 则其绝对误差：
1.1验并列为三大科学方法。
我们所学习的内容属于一门新学科——科学计算。
即现代意义下的计算数学。 主要研究在计算机上可计算的有效算法及相关理论。
本课程主要研究用计算机求解各种数学问题
的数值计算方法及其理论与软件实现
主要内容包括： 数值代数
Ax b
到x的前几位近似值a， 例如
x π 3.14159265
a1 0.00159265 取5位： a2 3.1416, a2 0.00000735
取3位： a1 3.14, 那么，它们的误差界的取法应为：
1 π 3.14 10 2 , 2
近似的，因此在此过程中，误差无处不在。
计 算 机 科 学 计 算 的 流 程 图
实际问题 数学模型 数值计算方法
编程实现算法 计算机数值结果
模型误差

方法误差或称为 截断误差 观测误差
舍入误差
1．模型误差 由实际问题抽象出数学模型，要简化许多 条件，这就不可避免地要产生误差．实际问题的解与数学模 型的解之间的误差
大连理工大学-矩阵与数值分析第一章上

大连理工大学-矩阵与数值分析第一章上


–对C的类型系统改进和扩充(更安全)
–支持面向对象
C++保持与C兼容(快速普及) C++不是纯粹的面向对象的语言

33
1.2 程序的编译过程
34
1.3 C++的词法记号
关键字 各种常量 操作符 标识符 分隔符

35
1.4 C++程序的结构
#include <iostream.h> int main() { cout<<”this is the start of something wonderful!”; cout<<endl; cout<<”And now we can say even more!”; return 0; }
Sub1
Sub2
….
Subn
各子流程实现----函数化 Func1 Func2 …. Funcn
根据系统的流程组建软件，通过函数的调用实现
17
面向对象思想
问题域 (Domain) 以问题域中的事物为中心思考问题 Object1 Object2
….
Objectn
对象归类----抽象化 Class1 Class2 …. Classn
返回类型
{ 函数体； }
50
函数名（形式参数1， 形式参数2，。。。，形式参数
3.2 参数的传递
值调用
#include <iostream.h> double Volume(double radius,double height); int main() { double v; v=Volume(3.0,3.0); cout<<"Volume="<<v<<endl; return 0; } double Volume(double radius,double height) { double result=3.14 * radius * radius * height; return result; }