共享知识分享快乐
大连理工大学
矩阵与数值分析上机作业
课程名称：矩阵与数值分析
研究生姓名：
12 交作业日时间：日20 月年2016
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
第1题
1.1程序：
Clear ;all n=input('请输入向量的长度n：') for i=1:n;
v(i)=1/i;
end
Y1=norm(v,1)
Y2=norm(v,2)
Y3=norm(v,inf)
1.2结果
n=10 Y1 =2.9290
Y2 =1.2449
Y3 =1
n=100 Y1 =5.1874
Y2 =1.2787
Y3 =1
n=1000 Y1 =7.4855
Y2 =1.2822
Y3 =1
N=10000 Y1 =9.7876
Y2 =1.2825
Y3 =1
1.3 分析
一范数逐渐递增，随着n的增加，范数的增加速度减小；二范数随着n的增加，逐渐趋于定值，无群范数都是1.
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
第2题
2.1程序
;clear all x(1)=-10^-15;dx=10^-18;L=2*10^3; i=1:L for
y1(i)=log(1+x(i))/x(i); d=1+x(i); d == 1if
y2(i)=1;else
y2(i)=log(d)/(d-1);end
x(i+1)=x(i)+dx;end x=x(1:length(x)-1););'r'plot(x,y1,on hold
plot(x,y2);
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
2.2 结果
2.3 分析
红色的曲线代表未考虑题中算法时的情况，如果考虑题中的算法则数值大小始终为1，这主要是由于大数加小数的原因。
第3题
3.1 程序
;clear all A=[1 -18 144 -672 2016 -4032 5376 -4608 2304 -512];x=1.95:0.005:2.05; i=1:length(x);for y1(i)=f(A,x(i)); y2(i)=(x(i)-2)^9;end figure(3);
plot(x,y1);;on hold
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
);'r'plot(x,y2,
F.m文件
y=f(A,x)function y=A(1); i=2:length(A);for y=x*y+A(i);;end3.2 结果
第4题
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
4.1 程序
;clear all n=input('请输入向量的长度n：')
A=2*eye(n)-tril(ones(n,n),0); i=1:n for A(i,n)=1;end n=length(A);U=A; e=eye(n);for i=1:n-1
[max_data,max_index]=max(abs(U(i:n,i))); e0=eye(n);
max_index=max_index+i-1; U=e0*U; e1=eye(n); j=i+1:n for
e1(j,i)=-U(j,i)/U(i,i);end
U=e1*U;中把变换矩阵存到P P{i}=e0;% L{i}=e1; e=e1*e0*e;end
k=1:n-2for Ldot{k}=L{k}; i=k+1:n-1for
Ldot{k}=P{i}*Ldot{k}*P{i};end
end Ldot{n-1}=L{n-1};LL=eye(n);PP=eye(n); i=1:n-1for PP=P{i}*PP;
LL=Ldot{i}*LL;end
b=ones(n,2);解方程 %b=e*b;x=zeros(n,1);x(n)=b(n)/U(n,n); i=n-1:-1:1for卑微如蝼蚁、坚强似大象．
共享知识分享快乐
x(i)=(b(i)-U(i,:)*x)/U(i,i);end计算逆矩阵%X=U^-1*e^-1*eye(n);AN=X'; result2{n-4,1}=AN;result1{n-4,1}=x;,n)'%d:\n'fprintf(fprintf('%d ',AN);
4.2 结果
n=5
1.0625 -0.875 -0.75 -0.5 -0.0625
-0.0625 0.0625 -0.75 1.125 -0.5
-0.0625 0.125 0.0625 1.25 -0.5
-0.0625 0.125
0.25 0.0625
1.5
0.0625
-0.5
-0.25
-0.0625 -0.125
n=10
1.0625 -0.875 -0.75 -0.5 -0.0625 1.0625 -0.875 -0.75 -0.5 -0.0625 -0.0625 1.125 0.0625 -0.75 -0.5 -0.5 0.0625 1.125 -0.75 -0.0625 -0.0625 0.0625 0.125 1.25 1.25 -0.0625 -0.5 0.0625 0.125 -0.5
-0.0625 0.25
0.25
0.0625 0.125
1.5 1.5 -0.0625 0.125
0.0625
0.0625 -0.0625 -0.125 -0.25 0.0625 -0.5 -0.0625 -0.125 -0.25 -0.5 -0.0625 -0.75 1.0625 -0.5 -0.0625 -0.875 -0.5 -0.75 1.0625 -0.875 -0.0625 -0.5 0.0625 1.125 -0.5 0.0625 1.125 -0.75 -0.0625 -0.75 1.25 0.125 0.0625 -0.0625 -0.0625 -0.5 -0.5 0.0625 0.125 1.25
0.25
-0.0625 -0.0625 1.5
0.125
0.0625 0.0625 0.25
0.125
1.5
-0.0625 -0.125 -0.25 0.0625
-0.5 0.0625 -0.0625 -0.125 -0.25
-0.5
同样的方法可以算出n=20,n=30时的结果，这里就不罗列了。
第5题
5.1 程序
;clear all n=input('请输入向量的长度n：10至20')
卑微如蝼蚁、坚强似大象．
共享知识分享快乐
i=1:n for j=1:n for
a(i,j)=1/(i+j-1);end
end j=1:n for sum=0; k=1:j-1for
sum = sum + l(j,k)^2;end
函数定义： % 目标函数 function f = fun(x) fm=0; for i=1:499 fmi = (1-x(i))^2 + 100*(x(i+1)-x(i)^2)^2; fm=fm+fmi; end f =fm; end % 梯度 function g = grad(x) g = zeros(500,1); g(1)=2*(x(1)-1)+400*x(1)*(x(1)^2-x(2)); for i=2:499 g(i)=2*(x(i)-1)+400*x(i)*(x(i)^2-x(i+1))+200*(x(i)-x(i-1)^2); end g(500) = 200*(x(500)-x(499)^2); end % 二阶梯度
function g = grad2(x) g = zeros(500,500); g(1,1)=2+400*(3*x(1)^2-x(2)); g(1,2)=-400*x(1); for i=3:500 g(1,i)=0; end for i=1:498 g(500,i)=0; end g(500,499)=-400*x(499); g(500,500)=200; for i=2:499 for j=1:500 if j==i-1 g(i,j)= -400*x(i-1); elseif j==i g(i,j)= 2+400*(3*x(i)^2-x(i+1))+200; elseif j==i+1 g(i,j)= -400*x(i); else g(i,j)=0; end end end end 1.最速下降法 function x_star = steepest(x0,eps) gk = grad(x0); res = norm(gk); k = 0; while res > eps && k<=50000 dk = -gk;
昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰
P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果：
求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1
数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )
则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
数值分析上机报告 院系： 学号： 姓名：
目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12)
作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比
Matlab作业3（数值分析） 机电工程学院（院、系）专业班组 学号姓名实验日期教师评定 1.计算多项式乘法(x2+2x+2)(x2+5x+4)。 答： 2. （1）将(x-6)(x-3)(x-8)展开为系数多项式的形式。（2）求解在x=8时多项 式(x-1)(x-2) (x-3)(x-4)的值。 答：（1） （2）
3. y=sin(x)，x从0到2π，?x=0.02π，求y的最大值、最小值、均值和标准差。 4.设ｘ＝[0.00.30.8 1.1 1.6 2.3]'，y=[0.500.82 1.14 1.25 1.35 1.40]'，试求二次多项式拟合系数，并据此计算x1=[0.9 1.2]时对应的y1。解：x=[0.0 0.3 0.8 1.1 1.6 2.3]'; %输入变量数据x y=[0.50 0.82 1.14 1.25 1.35 1.40]'; %输入变量数据y p=polyfit(x,y,2) %对x,y用二次多项式拟合,得到系数p x1=[0.9 1.2]; %输入点x1 y1=polyval(p,x1) %估计x1处对应的y1 p = -0.2387 0.9191 0.5318 y1 = a) 1.2909
5.实验数据处理：已知某压力传感器的测试数据如下表 p为压力值，u为电压值，试用多项式 d cp bp ap p u+ + + =2 3 ) ( 来拟 合其特性函数，求出a,b,c,d，并把拟合曲线和各个测试数据点画在同一幅图上。解： >> p=[0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9]; u=[10,11,13,14,17,18,22,24,29,34,39]; x=polyfit(p,u,3) %得多项式系数 t=linspace(0,10,100); y=polyval(x,t); %求多项式得值 plot(p,u,'*',t,y,'r') %画拟和曲线 x = 0.0195 -0.0412 1.4469 9.8267
大连理工大学《工程抗震》大作业
题目1：底部剪力法。 钢筋混凝土5层框架经质量集中后计算简图如下图所示，各层高均为3m ， 集中于各楼层的重力荷载代表值分别为： 1500kN G =，2550kN G =，3580kN G =，4600kN G =，5450kN G =。结构阻尼比0.05ξ=，自振周期为10.55s T =，Ⅰ1类 场地类别，设计地震分组为第一组，抗震设防烈度为8度（设计基本地震加速度为0.30g ）。按底部剪力法计算结构在多遇地震时的水平地震作用及地震剪力。 3580kN =2550kN =1500kN =（a ）计算简图 4600kN =5450kN = 解：查《建筑设计抗震规范》表5.1.4知，8度多遇地震，αmax=设计地震分组为第一组， Ι类场地，取Tg= Tg=＜T1=＜5Tg= α1=（Tg/T1）r η2αmax =（）××=≈ 查《建筑设计抗震规范》表5.2.1知，T 1=＞=×= 取δn=T1+=×+= 总水平地震作用标准值： F EK =α1Geq=×（500+550+580+600+450）×85%=
各楼层水平地震作用标准值： Fi=G i H i F EK (1-δn)/∑G j H j (i=1,2,3n) ∑G j H j =500×3 +550×6+580×9+600×12+450×15=23970KN ·m F 1=[500×3××]/23970= F 2=[550×6××]/23970= F 3=[580×9××]/23970= F 4=[600×12××]/23970= F 5=[450×15××]/23970= 计算各楼层的层间地震剪力 V 1= F 1+ F 2+ F 3+ F 4+ F 5=++++= V 2= F 2+ F 3+ F 4+ F 5=+++=152KN V 3= F 3+ F 4+ F 5=++= V 4= F 4+ F 5=+= V 5=F 5= 题目3：怎样判断土的液化如何确定土的液化严重程度，并简述抗液化措施。 答：饱和松散的砂土或粉土（不含黄土），地震时易发生液化现象，使地基承载力丧失或减弱，甚至喷水冒砂，这种现象一般称为砂土液化或地基土液化。其产生的机理为：地下水位以下的饱和砂土和粉土颗粒在地震作用下，土颗粒之间有变密的趋势。因空隙水不能及时排出，土颗粒就处于悬浮状态，形成如同液体一样的现象，即所谓的土的液化现象。地基土液化判别过程可以分为初步判断和标准贯入试验判别两大步骤。下面分别予以介绍。 1、初步判断 饱和的砂土或粉土（不含黄土）当符合下列条件之一时，可初步判别为不液化或不考虑液化影响： （1）地质年代为第四纪晚更新世（Q3）及其以前时且处于烈度7度或者8度地区时可判为不液化土。 （2）粉土的粘粒（粒径<0.005mm ）含量百分率当烈度为7度时大于10％、当烈度为8度时大于13％、当烈度为9度时大于16％，可判为不液化土。 （3）浅埋天然地基，当地下水位深度和覆盖非液化土层厚度满足下式之一时，可不考虑液化影响。 03w b d d d >+- 02 u b d d d >+-
2015.1.9 上机作业题报告 JONMMX 2000
1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果
1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序
高等数值分析上机作业
目录 上机作业1 (1) 上机作业2 (5) 上机作业3 (10) 上机作业4 (13) 上机作业5 (16) 上机作业6 (19) 上机作业7 (20) 上机作业8 (29)
第8章 函数逼近与曲线拟合 上机作业1： 最佳平方逼近 8-11.设()[]1,1,arcsin -∈=x x x f ， （1） 在{}32,,,1x x x span =φ中求()x f 的最佳平方逼近多项式； （2） 在{})(),(),(),(3210x T x T x T x T span =φ中求()x f 的最佳平方逼近多项式。 解：（1） 基于幂函数的最佳平方逼近 简单原理： 对于],[)(b a C x f ∈及一个线性无关函数组的集合 {}],,[)(,),(),(10b a C x x x span n ?=???φ 若存在,φ∈*S 使得 ()dx x S x f x x S x f x S x f b a x S x S ?-=-=-∈∈* 2)(2 2)(2 2 )]()([min )()(min )()(ρφ φ ，则称()x S *是 ()x f 在子集[]b a ,?φ中的最佳平方逼近函数。 取(),,,1,0,n j x x j j ==?就有{}n x x x span ,,,,12 =φ。对于任意的()φ∈x S ，有()∑==n j j j x a x S 0，()x S 为次数n ≤的多项式。 令)(x f 在},,,1{32x x x span =φ中的最佳平方逼近函数为 φ∈+++=3 *2**1*0*3 2)(x a x a x a a x S 通过求解法方程 ???? ? ? ? ??=??????? ????????? ??),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),(),() ,(321032103323130 3322212023121110130201000????????????????????????????????????f f f f a a a a 其中.arcsin )(,)(,)(,)(,1)(332210x x f x x x x x x x =====????
1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2;
绪论部分： 2、简述环境问题的分类？（10分） 答：环境问题是多方面的，但大致可分为两类：原生环境问题和次生环境问题。由自然力引起的为原生环境问题，也称为第一环境问题。由于人类生产和生活引起生态系统破坏和环境污染，反过来又危及人类自身和生存和发展的现象，为次生环境问题，也叫第二环境问题。原生环境问题和次生环境问题很难截然分开，它们之间常常存在着某种程度的因果关系和相互作用。 4、什么是环境化学，学习环境化学有什么意义？（10分） 答：环境化学是一门研究有害化学物质在环境介质中的存在、化学特性、行为和效应及其控制的化学原理和方法的科学。 意义：用来掌握污染来源，消除和控制污染，确定环境保护决策，以及提供科学依据诸方面都起着重要的作用。 5、简述环境化学的分支学科。（10分） 答：主要包括6类。 ①环境分析化学：是研究化学品的形态、价态、结构、样品前处理和痕量分析的学科。 ②环境污染化学：大气、水体和土壤环境化学，元素循环的化学过程。 ③污染控制化学：主要研究与污染控制有关的化学机制及工艺技术中化学基础性问题。 ④污染生态化学：是研究化学污染物在生态系统中产生生态效应的化学过程的学科。 ⑤环境计算化学：主要利用有效的数学近似以及电脑程序计算分子的性质。 ⑥环境生物化学：是研究环境化学品对生命影响的学科。 第一章： 1、地球环境主要由哪些圈层构成？英文单词？各之间有什么联系？各有哪些性 质？（10分） 答：地球环境主要由大气圈（atmosphere）、水圈（hydrosphere）、土壤圈（pedosphere）、岩石圈（lithosphere）和生物圈（biosphere）构成。 联系：大气圈、水圈、土壤圈和生物圈共同组成了地球环境系统，每个圈层都离不开
数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？
第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err
p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。
数值分析上机第四次作业 实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组 (1) 123421003131020141100155x x x x -?????? ? ? ?--- ? ? ?= ? ? ?-- ? ? ?-???? ?? 迭代20次或满足()(1) 1110k k x x --∞-<时停止计算。 (2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1 迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。 （3）*考虑模型问题，方程为 222222(),(,)(0,1)(0,1)(0,)1,(1,),01(,0)1,(,1),01 xy y x u u x y e x y D x y u y u y e y u x u x e x ??+=+∈=???==≤≤==≤≤ 用正方形网格离散化，若取1/,10h N N ==，得到100n =的线性方程组，并用共轭梯度法（CG 法）求解，并对解作图。 实验要求：迭代初值可以取01(,1,...,)ij u i j N ==，计算到32||||10k r -≤停止．本 题有精确解(,)xy u x y e =，这里k u 表示以k ij u 为分量的向量， u 表示在相应点(,)i j 上取值作为分量的向量． 实验一： （1） 编制函数子程序CGmethod 。 function [x,k]=CGmethod(A,b) n=length(A);x=zeros(n,1);r=b-A*x;rho=r'*r; k=0; while rho>10^(-12) & k<20 k=k+1; if k==1 p=r; else beta=rho/rho1; p=r+beta*p; end
2013级工科硕士研究生 《矩阵与数值分析》课程数值实验报告 大连理工大学 Dalian University of Technology
一、设 6 2 2 10 1 N N j S j = = - ∑，分别编制从小到大和从大到小的顺序程序分别计算 100001000000 , S S 并指出两种方法计算结果的有效位数。 程序代码： 从小到大： function f=s(N); %定义函数s f=0; %初始值为0 for j=N:-1:3 %j从3到n循环（从小到大） ft=1000000/(j^2-1); %Sj f=f+ft; %SN end 从大到小： function f=s(N); %定义函数s f=0; %初始值为0 for j=N:-1:3 %j从3到n循环（从小到大） ft=1000000/(j^2-1); %Sj f=f+ft; %SN end 执行结果： 从小到大： s(10000) ans = 4.16566671666167e+05 s(1000000) ans =
4.166656666671731e+05 有效数字：16,16 从大到小： s(10000) ans = 4.165666716661668e+05 s(1000000) ans = 4.166656666671667e+05 有效数字：16,16 分析： 小数和大数相加时，按照从大到小的顺序和按照从小到大的顺序得出的结果不同，前者由 于舍入误差的影响而使结果不准确，所以应避免大数吃小数的现象。 二、解线性方程组 1．分别利用Jacobi 迭代法和Gauss-Seidel 迭代法求解线性方程组Ax b =，其中常向量为()21n -维随机生成的列向量，系数矩阵A 具有如下形式 1111 11 1122n n n n n n n n T I I I A I I T I --------+-?? ?- ?= ? - ? -+? ? ， 其中1 211112n T --?? ? - ?= ?- ? -? ? 为1n -阶矩阵，1n I -为1n -阶单位矩阵，迭代法计算停止的条件为：10 12 10k k x x -+-<，给出10,100,1000n =时的不同迭代步数． 程序代码：
目录 一、设计任务书及原始数据 (2) 二、根据已知条件计算传动件的作用力 (3) 2.1计算齿轮处转矩T、圆周力F t 、径向力F r及轴向力F a .. 3 2.2计算链轮作用在轴上的压力 (3) 2.3计算支座反力 (4) 三、初选轴的材料，确定材料的机械性能 (4) 四、进行轴的结构设计 (5) 4.1确定最小直径 (5) 4.2设计其余各轴段的直径和长度，且初选轴承型号 (5) 4.3选择连接形式与设计细部结构 (6) 五．轴的疲劳强度校核 (6) 5.1轴的受力图 (6) 5.2绘制弯矩图 (7) 5.3绘制转矩图 (8) 5.4确定危险截面 (9) 5.5计算当量应力，校核轴的疲劳强度 (9) 六、选择轴承型号，计算轴承寿命 (10)
6.1计算轴承所受支反力 (10) 6.2计算轴承寿命 (11) 七、键连接的计算 (11) 八、轴系部件的结构装配图 (12) 一、设计任务书及原始数据 题目二：二级展开式斜齿圆柱齿轮减速器输出轴组合结构设计
表1 设计方案及原始数据 二、根据已知条件计算传动件的作用力 2.1计算齿轮处转矩T、圆周力F t、径向力F r及轴向力F a 已知：轴输入功率P=4.3kW，转速n=130r/(min)。 (1)齿轮上的力 转矩计算公式：T=9.550×106P/n 将数据代入公式中，得：T=315885(N·mm) 圆周力计算公式： F t =2T/，==416(mm) (认为是法面模数) 将转矩T带入其中，得：F t =1519(N) 径向力计算公式：F r =F t ×tanα/cos，= 将圆周力F t 带入其中，得：F r =558(N) 轴向力计算公式：F a = F t ×tan 将圆周力F t 带入其中，得：F a =216(N) 2.2计算链轮作用在轴上的压力 链轮的分度园直径 链速v= 链的圆周力F= 链轮作用在轴上的压力
大连理工大学《高层建筑结构》大作业 学习中心：浙江海盐奥鹏学习中心[23] 姓名：董刘明 学号：141350403661
题目二：风荷载计算题 某市的市区，有较多高层房屋的密集建筑群，欲建一座40层高的钢筋混凝土框架—剪力墙结构房屋，外形和质量沿高度方向基本呈均匀分布。房屋总高度140m H =，迎风面的房屋长50m L =。房屋总宽度30m B =，基本风压为 200.65kN m ω=。求楼顶处的风荷载标准值。 （按照教材及课件内容计算） 答：(1)风荷载体型系数 H/B=140/30=4.67>4 L/B=50/30=1.67>1.5 根据《高规》3.2.5条第4款3的规定，该楼已符合H/B>4的条件，但不 满足L/B ≯1.5的条件。故采用《高规》附录A 第1条规定的μs ： μs=0.8+(0.48+0.03H/L)=0.8+0.48+0.03×140/50=1.364 (2)该房屋建于市区，因指明是属于房屋较高，建筑群密集的市区，根据《高规》3.2.3条的规定，地面粗糙度属于D 类。 Ζ=140m βz=1+(ρz §ν/μz)=1+Ζ§ν/H μz=1+[(140×1.491×0.487)/140×1.50]=1.484 作用于屋顶处的风荷载标准值： ωk=βzμsμzω0=1.484×1.364×1.50×0.65=1.974kN/m2 题目四：整体墙、小开口整体墙、联肢墙、带刚域框架和单独墙肢等计算方法 的特点及适用条件是什么？ 答：整体墙计算中认为平面假定仍然适用，截面中的正应力符合直线分布规律，按照整体悬臂墙方法计算墙在水平荷载作用下的截面内力和变形。适用于墙面上门窗、洞口等开孔面积不超过墙面面积的15‰，且洞口间净距及孔洞至墙边的净距大于孔洞边长尺寸。 小开口整体墙的整体弯曲变形仍是主要的，其内力可按材料力学公式计算，但需考虑局部弯曲的影响稍作修正。适用于洞口稍大一些，且洞口沿竖向成列布置，洞口面积超过墙面总面积的16％，但洞口对剪力墙的受力影响较小。 联肢墙计算比较复杂，在引入连杆连续化假定，对于竖向布置规则的双肢墙
数值分析上机报告 《 院系： 学号： 姓名：
目录 作业1、舍入误差与有效数 (1) 1、函数文件 (1) 2、函数文件 (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (2) 2、局部收敛性 (3) （1）最大δ值文件 (3) （2）验证局部收敛性 (4) 3、心得 (5) 作业3、列主元素Gauss消去法 (6) 1、列主元Gauss消去法的通用程序 (6) 2、解题中线性方程组 (7) 3、心得 (8) 作业4、三次样条插值函数 (8) 1、第一型三次样条插值函数通用程序： (8) 2、数据输入及计算结果 (10)
作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序11 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么 程序： 1、函数文件 function S=cxdd(N) S=0; i=; while (i<=N) S=S+(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 2、函数文件 function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 3、两种方法有效位数对比 精确值函数： function S=jqz(N) S=*运行结果（省略>>）
数值分析思考题1 1、讨论绝对误差（限）、相对误差（限）与有效数字之间的关系。 2、相对误差在什么情况下可以用下式代替？ 3、查阅何谓问题的“病态性”，并区分与“数值稳定性”的不同点。 4、 取 ，计算 ，下列方法中哪种最好？为什么？ （1）()3 322 -，（2）() 2 752-，（3） () 3 1 322+，（4） () 6 121 +，（5）99702- 数值实验 数值实验综述：线性代数方程组的解法是一切科学计算的基础与核心问题。求解方法大致可分为直接法和迭代法两大类。直接法——指在没有舍入误差的情况下经过有限次运算可求得方程组的精确解的方法，因此也称为精确法。当系数矩阵是方的、稠密的、无任何特殊结构的中小规模线性方程组时，Gauss 消去法是目前最基本和常用的方法。如若系数矩阵具有某种特殊形式，则为了尽可能地减少计算量与存储量，需采用其他专门的方法来求解。 Gauss 消去等同于矩阵的三角分解，但它存在潜在的不稳定性，故需要选主元素。对正定对称矩阵，采用平方根方法无需选主元。方程组的性态与方程组的条件数有关，对于病态的方程组必须采用特殊的方法进行求解。 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 r e x x e x x ***** -== 2141.≈()6 21-
大连理工大学优化方法上机 作业 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
优化方法上机大作业 学院：电子信息与电气工程学部 姓名： 学号： 指导老师：
上机大作业（一） %目标函数 function f=fun(x) f=100*(x(2)-x(1)^2)^2+(1-x(1))^2; end %目标函数梯度 function gf=gfun(x) gf=[-400*x(1)*(x(2)-x(1)^2)-2*(1-x(1));200*(x(2)-x(1)^2)]; End %目标函数Hess矩阵 function He=Hess(x) He=[1200*x(1)^2-400*x(2)+2,-400*x(1); -400*x(1), 200;]; end %线搜索步长 function mk=armijo(xk,dk) beta=0.5; sigma=0.2; m=0; maxm=20; while (m<=maxm) if(fun(xk+beta^m*dk)<=fun(xk)+sigma*beta^m*gfun(xk)'*dk) mk=m; break; end m=m+1; end alpha=beta^mk newxk=xk+alpha*dk fk=fun(xk) newfk=fun(newxk) %最速下降法 function [k,x,val]=grad(fun,gfun,x0,epsilon) %功能：梯度法求解无约束优化问题：minf(x) %输入：fun，gfun分别是目标函数及其梯度，x0是初始点， % epsilon为容许误差 %输出：k是迭代次数，x，val分别是近似最优点和最优值 maxk=5000; %最大迭代次数
中特作业 中特课学习体会 学院（系）：化工与环境生命学部专业：环境科学与工程 学生姓名：牛丹阳 学号：21607331 完成日期：2017年6月5日 大连理工大学 DalianUniversity of Technology
通过中特课程的学习，我既学习到了我党政治道路的有关理论知识，同时又学习到了如何将理论联系实际，将所学到的知识应用到国家治理的实际中去，走好中国特色的社会主义道路。此次作业，我将着重探讨我对于“一带一路”的学习体会。 一带一路是指“丝绸之路经济带”和“21世纪海上丝绸之路”，“一带一路”不是一个实体和机制，而是合作发展的理念和倡议，是充分依靠中国与有关国家既有的双多边机制，借助既有的、行之有效的区域合作平台，旨在借用古代“丝绸之路”的历史符号，高举和平发展的旗帜，积极主动地发展与沿线国家的经济合作伙伴关系，共同打造政治互信、经济融合、文化包容的利益共同体、命运共同体和责任共同体。 自一带一路战略实施以来，取得了丰硕的成果。亚投行、蒙内铁路、卫星通讯、中匈协议使得中国与丝路沿途国家得以分享优质产能，共商项目投资、共建基础设施、共享合作成果。这其中，老师在课上讲的中巴铁路给我留下了深刻的印象。 中巴铁路的修建是具有极为重要的历史意义的。主要是以下几点： 1、我们自己消化产能过剩，拓展发展空间，寻求陆地出口，获得印度洋出海口投资、商贸、军事等战略利益，加强中国在地缘政治中的地位。 2、中巴铁路连接出海口瓜达尔港，瓜达尔港有着极其特殊的军事价值。最近中国承建的瓜达尔港开航典礼刷爆了朋友圈，这是迄今为止，中国“一带一路”战略取得的最重要成果，也有人形容中国获得美苏没有抢到的地球支点。 中国一旦建设好了瓜达尔港这个深水良港，中国海军将控制中东石油运输的大动脉赫尔木兹海峡，这是印度的能源生命线，只要印度轻易妄动，中国和巴基斯坦随时可以切断印度的石油供应。再者，只要中国在瓜达尔港建类似于萨德那样的情报监测系统，那么整个中东和南亚都没有什么重要的目标能逃过中国天眼。