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 =时的不同迭代步数. 程序代码:
Jacobi迭代法
function [u,l,b,x,kk,delta,A]=Ja(n)
t=2.*eye(n-1);
s=-eye(n-2);
z=zeros(n-2,1); %生成n-1行1列零向量
z1=zeros(1,n-1);
s1=[z s];
s1=[s1;z1];
s2=[s z];
s2=[z1;s2];
T=t+s1+s2; %生成T
i=-eye((n-1)*(n-2));
a1=[zeros(((n-1)*(n-2)),n-1) i];
a2=[a1;zeros(n-1,((n-1)*(n-1)))];
b1=[i zeros(((n-1)*(n-2)),n-1)];
b2=[zeros(n-1,((n-1)*(n-1)));b1];
d=4.*eye((n-1)*(n-1)); %生成D
b=round(randn((n-1)*(n-1),1)); %生成b
A0=[];
for k=1:1:n-1
x=T+2.*eye(n-1);
A0=blkdiag(A0,x); %生成A
end
A=A0+a2+b2; %生成A
u=-(triu(A)-d); %生成U
l=-(tril(A)-d); %生成L
x=randn((n-1)*(n-1),1); %生成初始x
%%以下为计算迭代部分
B=inv(d)*(l+u);
ff=inv(d)*b;
kk=1;
derta=10;
delta=[];
x0=x;
while derta>10^-10&(kk<800);
x0=x;
x=B*x+ff;
derta=norm(x-x0);
kk=kk+1;
delta=[delta derta];
end
plot(delta);
End
执行结果:
运行结果:[u,l,b,x,kk,delta,A]=Ja(10)情况
图1 迭代曲线收敛步数465 收敛值9.899081611039407e-11
[u,l,b,x,kk,delta,A]=Ja(100)情况下
图2 迭代步数48988 收敛终值9.997187164447103e-11 当n=1000时,矩阵规模超出内存。
程序代码:
Gauss-Seidel迭代法
function kk=Gs(n)
t=2.*eye(n-1);
s=-eye(n-2);
z=zeros(n-2,1); %生成n-1行1列零向量
z1=zeros(1,n-1);
s1=[z s];
s1=[s1;z1];
s2=[s z];
s2=[z1;s2];
T=t+s1+s2; %生成T
i=-eye((n-1)*(n-2));
a1=[zeros(((n-1)*(n-2)),n-1) i];
a2=[a1;zeros(n-1,((n-1)*(n-1)))];
b1=[i zeros(((n-1)*(n-2)),n-1)];
b2=[zeros(n-1,((n-1)*(n-1)));b1];
d=4.*eye((n-1)*(n-1)); %生成D
b=round(randn((n-1)*(n-1),1)); %生成b
A0=[];
for k=1:1:n-1
x=T+2.*eye(n-1);
A0=blkdiag(A0,x); %生成A
end
A=A0+a2+b2; %生成A
u=-(triu(A)-d); %生成U
l=-(tril(A)-d); %生成L
x=randn((n-1)*(n-1),1); %生成初始x
%%以下为Gauss迭代法部分
B=(d-l)\u;
F=(d-l)\b;
x0=x;
derta=10;
kk=0;
delta=[];
while derta>10^-10&(kk<800);
x0=x;
x=B*x+F;
derta=norm(x-x0);
kk=kk+1;
delta=[delta derta];
end
plot(delta);
End
执行结果:
f=Gs(10)
图3 G-S 法n=10 迭代步数220 收敛终值9.623878789391471e-11
f=Gs(100) 分析:
Gauss-Seidel 迭代法是Jacobi 迭代法的改进,求解该线性方程组,Gauss-Seidel 迭代法和Jacobi 迭代法都收敛,并且Gauss-Seidel 迭代法收敛速度比Jacobi 迭代法快。但是具体的问题要具体分析。
2. 用Gauss 列主元消去法、QR 方法求解如下方程组:
123412110250
34.
179212812
18x x x x -?????? ?
? ? ? ? ?= ? ? ? ?
? ?---????
??
程序代码:
Gauss列主元消去法
function x=Gsl(A,b)
%%输入系数矩阵A和常数列向量b,高斯列主元法
Ab=[A b]; %构成增广矩阵
[m,n]=size(Ab);
for i=1:1:m
g=tril(Ab);
g=abs(g);
mx=find(g(:,i)==max(g(:,i))); %找到第m列最大值所在行数为mx
Ab([i mx],:)=Ab([mx i],:); %将第m行与mx行交换for k=(i+1):1:m
xb=-Ab(k,i)/Ab(i,i);
xc=xb*Ab(i,:);
Ab(k,:)=Ab(k,:)+xc;
end
end
%%至此高斯列主元结束,下面为回带法
x=zeros(m,1);%初始化解向量x
x=x+1;
A=Ab(:,1:m);
b=Ab(:,n);
for i=m:-1:1%回带法求x向量
x(i,1)=0;
format long; %取长精度
x(i,1)=(b(i,1)-A(i,:)*x)./A(i,i);
end
End
执行结果:
x=(-1 0 1 2)T
程序代码:
QR 法
function [Q,R,x,y]=qrf(A,b) [Q,R]=qr(A);%求A 的QR 分解矩阵 y=Q\b;%利用回带法求解 x=R\y; End 执行结果:
x=(-1 0 1 2)T
y = (5.2590 -13.6412 -3.1085 0.7734)T
分析:
从上边的编程可以看出Gauss 列主元消去法与QR 方法求解的结果相同。
三、非线性方程的迭代解法
1.求方程
()222sin ln 160x f x e x x x =++--=
的根,迭代停止的条件为:10110k k x x -+-<;
程序代码:
%% 二分法迭代求方程的根,迭代区间选为[1,3] clear all; syms x;
f=exp(x)+2*x^2+2*sin(x)-log(x)-16; %函数f(x)的表达式 i=0; %二分次数记数 a=1; %求根区间左端 b=3; %求根区间右端 e=1e-10;
fa=subs(f,x,a); fb=subs(f,x,b);
c=(a+b)/2; %计算区间中点
fc=subs(f,x,c);
while abs(fc)>=e; %判断f(c)是否为零点
if fa*fc>=0; %判断左侧区间是否有根
fa=fc;
a=c;
else
fb=fc;
b=c;
end
c=(a+b)/2;
fc=subs(f,x,c);
i=i+1;
end
fprintf('采用二分法迭代求解方程的根为:%.10f \n',c);
fprintf('采用二分法迭代求解方程的根的迭代次数:%d \n',i); 执行结果:
采用二分法迭代求解方程的根为:1.9629226832
采用二分法迭代求解方程的根的迭代次数:35
2.利用Newton迭代法求多项式
432
--+-
x x x x
360
311=
的所有实零点,注意重根的问题。
程序代码:
%% Newton迭代法求多项式的所有实根
clear all;
syms x;
f=x^4-3*x^3-3*x^2+11*x-6; % 多项式f(x)
fx=diff(f); %对fx求导
fxx=diff(fx);
%x0=-1.5; %选取初始点
x0=0.5;
%x0=2.5;
f0=subs(f,x,x0); fx0=subs(fx,x,x0); e=1e-5; k=0;
y=x0-f0/fx0;
while abs(y-x0)>=e x0=y;
f0=subs(f,x,x0); fx0=subs(fx,x,x0); y=x0-f0/fx0; k=k+1; end
if (subs(fxx,x,y)~=0)&&(subs(fx,x,y)<=1e-4)
fprintf('用Newton 迭代法所求得的多项式的根是2重根:\n'); fprintf('根为:%.4f \n,',y); fprintf('迭代次数为:%d \n',k); else
fprintf('用Newton 迭代法所求得的多项式的根是单根:\n'); fprintf('根为:%.4f \n,',y); fprintf('迭代次数为:%d \n',k); end
执行结果:
(1).用Newton 迭代法所求得的多项式的根是2重根:根为:1.0000 ,迭代次数为:15 (2). 用Newton 迭代法所求得的多项式的根是单根:根为:-2.0000 ,迭代次数为:5 (3). 用Newton 迭代法所求得的多项式的根是单根:根为:3.0000 ,迭代次数为:7
分析:
本程序利用牛顿法迭代法进行迭代计算,公式为)
()
('
1k k k k x f x f m
x x -=+, 另外牛顿迭代法的初始点应该在零点附近,否则可能迭代不收敛。此外因为在x =1是二重根,需要判断下在x=1点处的导数。
四、数值积分
分别用三点Gauss 型求积公式计算积分
1
x e dx -?
程序代码:
clear all; syms x; f=exp(-x);
L=0.5*(5*x^3-3*x); % 三次Legendre 多项式 a=solve(L); b=a'; for i=1:3
A(1,i)=2/((1-b(1,i)^2)*[subs(diff(L),x,b(1,i))]^2); end I=0; for i=1:3
I=I+0.5*A(1,i)*subs(f,x,0.5*b(1,i)+0.5); end
fprintf(‘采用Gauss 型求积公式求得的结果为:’) disp(double(I));
执行结果:
采用Gauss 型求积公式求得的结果为: 0.6321
五、插值与逼近
1.给定[]1,1-上的函数()2
1
1f x x =
+,请做如下的插值逼近: (1)构造等距节点分别取5=n ,8=n ,10=n 的Lagrange 插值多项式; (2)取Chebyshev 多项式()()x n x T n arccos cos ?=的零点:
πn
k x k 21
2cos
-=,n k ,,1 =
作插值节点构造10=n 的插值多项式
()x f 和上述的插值多项式均要求画出曲线图形(用不同的线型或颜色表示不同的曲线)。
(1) 程序代码:
clear all;
n=5;%不同的n 代表不同的插值点数,可以修改 for i=1:n+1
xdata(i)=-1+(2/n)*(i-1); fdata(i)=1/(1+xdata(i)^2); end
f=Lagrange(xdata,fdata);%Lagrange 插值函数 tmp=100; for i=1:tmp+1
xxdata(i)=-1+(2/tmp)*(i-1); ffdata(i)=1/(1+xxdata(i)^2); end
ff=subs(f,xxdata); figure(1);
plot(xxdata,ffdata); hold on; grid on;
plot(xxdata,ff,'r');
title('红色且在-0.8附近值较小的曲线表示的是插值多项式,蓝色的是原函数f(x)'); xlabel('n=5时的曲线');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Lagrange 插值多项式求解
%xdata:插值点,fdata:插值点的函数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function f=Lagrange(xdata,fdata) syms x; f=0;
npl=length(xdata);
for k=1:1:npl
index=setdiff(1:npl,k);
f=f+fdata(k)*prod((x-xdata(index))./(xdata(k)-xdata(index)));
end
执行结果:
仿真曲线:
在-0.8附近值较小的曲线表示的是插值多项式,另一个是原函数f(x)'的曲线n=5
n=8
n=10
分析:
从编程可以得出插值的点数越多,得到的插值函数与原函数的误差越小。
(2)
程序代码:
clc;
clear all;
n=10;%不同的n代表不同的插值点数,可以修改
for i=1:1:n
xdata(i)=cos(((2*i-1)/(2*n))*pi);
fdata(i)=1/(1+xdata(i)^2);
end
f=Lagrange(xdata,fdata);%Lagrange插值函数
tmp=100;
for i=1:1:tmp
xxdata(i)=cos(((2*i-1)/(2*tmp))*pi);
ffdata(i)=1/(1+xxdata(i)^2);
end
ff=subs(f,xxdata);
figure(1);
plot(xxdata,ffdata);
hold on;
grid on;
plot(xxdata,ff,'r*');
title('红色且用"*"表示的曲线是插值多项式,蓝色的是原函数f(x)');
xlabel('n=10时的曲线');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %Lagrange插值多项式求解
%xdata:插值点,fdata:插值点的函数值
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function f=Lagrange(xdata,fdata)
syms x;
f=0;
npl=length(xdata);
for k=1:1:npl
index=setdiff(1:npl,k);
f=f+fdata(k)*prod((x-xdata(index))./(xdata(k)-xdata(index)));
End
执行结果:
仿真的多项式曲线:
3.观察物体的直线运动,得到如下数据
求运动学方程s at bt
=+。
程序代码:
%% 最小二乘拟和二次曲线
clear all;
syms x;
f=[1;x;x^2];
t=[0,0.9,1.9,3.0,3.9,5.0]';
s=[0,10,30,51,80,111]';
A=zeros(3,3);
b=zeros(3,1);
%% 构造系数矩阵A
for i=1:3
for j=1:3
sum_A=0;
for k=1:6
sum_A=sum_A+subs(f(i,1),x,t(k))*subs(f(j,1),x,t(k));
k=k+1;
end
A(i,j)=sum_A;
end
end
%% 构造b
for i=1:3
sum_b=0;
for k=1:6
sum_b=sum_b+subs(f(i,1),x,t(k))*s(k);
k=k+1;
end
b(i,1)=sum_b;
end
%% 待定系数求解
a=inv(A)*b;
disp('采用最小二乘拟合二次曲线的系数分别为:')
disp(a);
执行结果:
采用最小二乘拟合二次曲线的系数分别为:
-0.6179
11.1601
2.2684
分析:
从表达式可以看出,当每个点带入求出的公式所求的值与给出的点的值不相等,也就是有一定的误差,但是最小二乘法能够尽可能地表示运动的方程的趋势,即均方误差最小。
大连理工大学山上礼堂常用数据一览舞台 舞台上方横幅尺寸:14M*1M, 在12M范围内刻字(相应的舞台宽是14.2M) 舞台两侧台口竖条长:7.2M。 舞台背景喷绘尺寸:13M*6.5M (12M*6M) 后台左右两扇门的尺寸:130*225 cm 后台两侧的横梁:2.9M 舞台上左右两个音箱的尺寸:115*60 cm 从观众席向背景喷绘方向,左右两边依次是 红幕 绿幕1 绿幕2 绿幕3 粉幕 背景喷绘 其中:绿幕1紧贴红幕;绿幕3紧贴粉幕 红幕——圆弧形舞台边缘的最远点:3M 绿幕1——绿幕2距离: 175cm 绿幕2——粉幕距离:240cm 绿幕3——背景喷绘距离:680cm 观众席 观众席2楼(舞台对面)横幅15M 观众席两边竖条幅(即XX学院祝大会圆满成功的位置)尺寸:0.9*7.5M 一楼观众席,俯视的话可以分成六个区域 舞台
123 456 调音台 区域一:14排12列161个座位 区域二:14排17列251个座位(678排嘉宾席49个座)区域三:14排12列161个座位 区域四:10排12列120个座位 区域五:10排17列165个座位 区域六:10排12列120个座位 二楼观众席,俯视的话可以分成四个区域 舞台 12 34 调音台 区域一:6排22列115个座位 区域二:6排22列114个座位 区域三、四:8排22列400个座位 注:区域边缘呈锯齿状 前厅 前厅两侧宣传栏尺寸1.14M*3.94M 前厅柱子间距4.93M 前厅瓷砖壁画尺寸6.2M * 2.4M 礼堂正门 注:礼堂正面有四个竖直的突出部分,称为“柱子” 楼前中间柱子之间的间距5.8M 楼前两边柱子之间的间距12.4M
第一题 Lagrange插值函数 function y=lagrange(x0,y0,x); n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end x0=[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+12 x=[1:0.1:10]; x=x'; plot(x0,y0,'r'); hold on plot(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; end d1=sqrt(d1); d2=sqrt(d2); d3=sqrt(d3); fun=[f1 f2 f3;d2 d3 d4]; return;
大连理工大学2017年研究生矩阵与数值分析考试 考试日期:2017年6月5日 一、填空题(50分,每空2分) 1.a=0.3000经过四舍五入具有4位有效数字,则 x a a -≤,ln ln x a -≤ 2.已知X=(1,5,12)T ,Y=(1,0,a)T ,则由X 映射到Y 的Householder 矩阵为:,计算||H||2=,cond 2(H)= 3.根据3次样条函数的性质(后面-前面=a (x-x0)3),一个求其中的参数b== 4.2 '3u u t =,写出隐式Euler 格式: 梯形法格式: 5.已知A=XX T ,其中X 为n 维列向量,则||A||2=,||A||F =,矩阵序列的极限:2lim k k A A →∞?? ? ? ?? = 6.A=LU ,其解为x ,写出一步迭代后的改善格式: 7. 531A -?? ? = ? ?-?? ,请问通过幂法与反幂法计算出的特征值分别是, 8.1111A ?? ?= ? ??? ,sin A =,823A A A +-=,At e =,d d At e t =,2 1At e dt ?= 9. ()()()()2 1 2 012f x dx A f A f A f =++?是Newton-cotes 公式,则1 A =,具有代数精度= 10. f(x)=7x 7+6x 6+…+x ,f[20,21,22….,28]= 11. 0.40.200.5A ??= ???,1 k k A ∞=∑= 12.f(0)=1,f(1)=-1,f(2)=1,f(3)=19,请问对该节点进行插值后最高次的系数= 还有2空没有回忆出来,但是比上面题目还简单,因此不用担心。 二、121232352A -?? ?=-- ? ?--??,121b ?? ? = ? ?-?? (1)计算LU 分解 (2)利用LU 求逆矩阵 (3)写出G-S 格式(12分)
大连理工大学《工程抗震》大作业
题目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 >+-