矩阵与数值分析.docx
- 文档编号:5482765
- 上传时间:2023-05-08
- 格式:DOCX
- 页数:25
- 大小:147.53KB
矩阵与数值分析.docx
《矩阵与数值分析.docx》由会员分享,可在线阅读,更多相关《矩阵与数值分析.docx(25页珍藏版)》请在冰点文库上搜索。
矩阵与数值分析
2013级工科硕士研究生
《矩阵与数值分析》课程数值实验题目
一、设
,分别编制从小到大和从大到小的顺序程序分别计算
并指出两种方法计算结果的有效位数。
Matlab程序如下:
function[si,sd]=S(N)
formatlong;
si=0;sd=0;
forj=N:
-1:
2
si=1.0e6/(j^2-1)+si;
end
forj=2:
N
sd=1.0e6/(j^2-1)+sd;
end
end
在matlab命令窗口中输入:
[si,sd]=S(10000)
运行结果:
si=7.499000049995000e+005
sd=7.499000049994994e+005
在matlab命令窗口中输入:
[si,sd]=S(1000000)
运行结果:
si=7.499990000005000e+005
sd=7.499990000005200e+0051
结果分析:
si为从大到小的顺序求和的值,sd为从小到大的顺序求和的值。
当N分别为10000和1000000时,si分别为7.499000049995000e+005和7.499990000005000e+005,可以看出这两个数的有效值均为13位;而sd分别为7.499000049994994e+005和7.499990000005200e+005,这两个数的有效值均为16位。
这就出现了我们在矩阵理论课上所学的“大数吃小数”的问题。
为了使结果更为精确我们必须避免在四则运算中出现“大数吃小数”的情况,应该按从小到大的顺序进行求和。
二、解线性方程组
1.分别利用Jacobi迭代法和Gauss-Seidel迭代法求解线性方程组
,其中常向量为
维随机生成的列向量,系数矩阵
具有如下形式
,
其中
为
阶矩阵,
为
阶单位矩阵,迭代法计算停止的条件为:
,给出
时的不同迭代步数.
求解系数矩阵A和向量b的Matlab程序如下:
function[A,b,x0]=jz(n)
fori=1:
n-1%此处n赋值即可,如n=100
forj=1:
n-1
if(i==j)
T(i,j)=2;
end
if(i==(j+1))
T(i,j)=-1;
end
if(i==(j-1))
T(i,j)=-1;
end
end
end
I=eye(n-1);
k=1;
formm=1:
(n-1)
A(k:
(k+n-2),k:
(k+n-2))=T+2*I;
k=k+n-1;
end
k=1;
forxx=1:
(n-2)
A(k:
(k+n-2),(k+n-1):
(k+2*n-3))=-I;
A((k+n-1):
(k+2*n-3),k:
(k+n-2))=-I;
k=k+n-1;
end
b=randn((n-1)^2,1);
x0=zeros((n-1)^2,1);
Jacobi迭代法Matlab程序如下:
functionn=jacobi(A,b,x0)
eps=1.0e-10;
M=10000;
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
B=D\(L+U);
f=D\b;
x=B*x0+f;
n=1;%迭代次数
whilenorm(x-x0)>=eps
x0=x;
x=B*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛!
');
return;
end
end
Gauss-Seidel迭代法Matlab程序如下:
functionn=gauseidel(A,b,x0)
eps=1.0e-10;
M=10000;
D=diag(diag(A));%求A的对角矩阵
L=-tril(A,-1);%求A的下三角阵
U=-triu(A,1);%求A的上三角阵
G=(D-L)\U;
f=(D-L)\b;
x=G*x0+f;
n=1;%迭代次数
whilenorm(x-x0)>=eps
x0=x;
x=G*x0+f;
n=n+1;
if(n>=M)
disp('Warning:
迭代次数太多,可能不收敛!
');
return;
end
end
在matlab命令窗口中输入:
[A,b,x0]=jz(10);
J10=jacobi(A,b,x0)
G10=gauseidel(A,b,x0)
[A,b,x0]=jz(20);
J20=jacobi(A,b,x0)
G20=gauseidel(A,b,x0)
[A,b,x0]=jz(30);
J30=jacobi(A,b,x0)
G30=gauseidel(A,b,x0)
运行结果:
J10=436;G10=214
J20=1810;G20=913
J30=3990;G30=2001
结果分析:
N=10且M=1000时,Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为436和214;N=20且M=1000时,Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为1810和913;N=30且M=1000时,Jacobi和Gauss-seidel算法的迭代次数分别为3990和2001次。
从以上结果可知在该问题下Jacobi算法的收敛性没有Gauss-seidel算法的收敛性好,原因在于Jacobi算法迭代过程中,对已算出的信息未加充分利用,一般来说,后面计算的值要比前面的计算值要精确些,而Gauss-seidel算法则充分利用了已求出来的信息,所以此算法的收敛性更为好一些。
2.用Gauss列主元消去法、QR方法求解如下方程组:
Gauss列主元消去法Matlab程序如下:
function[x,XA]=GaussXQLineMain(A,b)
N=size(A);
n=N
(1);
index=0;
fori=1:
(n-1)
me=max(abs(A(1:
n,i)));%选取列主元
fork=i:
n
if(abs(A(k,i))==me)
index=k;%保存列主元所在的行
break;
end
end
temp=A(i,1:
n);
A(i,1:
n)=A(index,1:
n);
A(index,1:
n)=temp;
bb=b(index);
b(index)=b(i);
b(i)=bb;%交换主行
forj=(i+1):
n
if(A(i,i)==0)
disp('对角元素为0!
');
return;
end
l=A(j,i);
m=A(i,i);
A(j,1:
n)=A(j,1:
n)-l*A(i,1:
n)/m;
b(j)=b(j)-l*b(i)/m;%消元
end
end
fori=n:
-1:
1
if(i s=A(i,(i+1): n)*x((i+1): n,1); else s=0; end x(i,1)=(b(i)-s)/A(i,i); end XA=A; QR法Matlab程序如下: functionx=qrxq(A,b) N=size(A); n=N (1); B=A;%保存系数矩阵 A(1: n,1)=A(1: n,1)/norm(A(1: n,1));%将A的第一列正规化 fori=2: n forj=1: (i-1) A(1: n,i)=A(1: n,i)-dot(A(1: n,i),A(1: n,j))*A(1: n,j); %使A的第i列与前面所有的列正交 end A(1: n,i)=A(1: n,i)/norm(A(1: n,i)); %将A的第i列正规化 end Q=A;%分解出来的正交矩阵 R=transpose(Q)*B; bb=transpose(Q)*b; fori=n: -1: 1 if(i s=R(i,(i+1): n)*x((i+1): n,1); else s=0; end x(i,1)=(bb(i)-s)/R(i,i); end 在matlab命令窗口中输入: A=[12-11;2503;1792;8-1-21]; b=[0;4;12;-8]; x_G=GaussXQLineMain(A,b) x_QR=qrxq(A,b) 运行结果: x_G=-1.00000000000000 -0.00000000000000 1.00000000000000 2.00000000000000 x_QR=-1.00000000000000 -0.00000000000001 1.00000000000000 2.00000000000002 三、非线性方程的迭代解法 1.求方程 的根,迭代停止的条件为: ; 使用弦截法求解此方程的根,弦截法的Matlab程序如下: functionroot=Secant(f,a,b) %f: 方程; %a: 区间左端点; %b: 区间右端点; %root: 求得的根; eps=1.0e-10; f1=subs(sym(f),findsym(sym(f)),a); f2=subs(sym(f),findsym(sym(f)),b); if(f1==0) root=a; end if(f2==0) root=b; end if(f1*f2>0) disp('两端点函数值乘积大于0! '); return; else tol=1; fa=subs(sym(f),findsym(sym(f)),a); fb=subs(sym(f),findsym(sym(f)),b); root=a-(b-a)*fa/(fb-fa); while(tol>eps) r1=root; fx=subs(sym(f),findsym(sym(f)),r1); s=fx*fa; if(s==0) root=r1; else if(s>0) root=b-(r1-b)*fb/(fx-fb); else root=a-(r1-a)*fa/(fx-fa); end end tol=abs(root-r1); end end 利用matlab的绘图功能可以确定求根区间为(1,3) 在matlab命令窗口中输入: root=Secant('exp(x)+2*x^2+2*sin(x)-log(x)-16',1,3) 运行结果: root=1.96292268316446 2.利用Newton迭代法求多项式 的所有实零点,注意重根的问题。 Newton迭代法的Matlab程序如下: x=-3;d=-1; whileabs(d-x)>0.5*10^(-8)%在区间(-3,-1) d=x; x=x-(x^4-3*x^3-3*x^2+11*x-6)/(4*x^3-9*x^2-6*x+11);%牛顿法 end x0=d x=0;d=2; whileabs(d-x)>0.5*10^(-8)%在区间(0,2)的重根 d=x; x=x-2*(x^4-3*x^3-3*x^2+11*x-6)/(4*x^3-9*x^2-6*x+11); end x1=d x=2.5;d=4; whileabs(d-x)>0.5*10^(-8)%在区间(2.5,4) d=x; x=x-(x^4-3*x^3-3*x^2+11*x-6)/(4*x^3-9*x^2-6*x+11); end x2=d 运行结果: x0=-2.00000000151233 x1=1.00000000065769 x2=3.00000000000002 结果分析: 由以上结果可知,newton迭代法所得结果与初始值的选取密切相关。 不同的初始值会产生不同的结果,甚至会影响牛顿迭代法收敛性。 四、数值积分 用三点Gauss型求积公式计算积分 n点Gauss型求积公式matlab程序如下: functionI=question4(f,a,b,n) %f为求积方程; %a,b分别为积分上下限; %n为Gauss点个数; [XK,AK]=grule(n); XK1=((b-a)/2)*XK+((a+b)/2); I=((b-a)/2)*sum(AK.*subs(f,findsym(f),XK1)); function[bp,wf]=grule(n) %[bp,wf]=grule(n) bp=zeros(n,1);wf=bp;iter=2;m=fix((n+1)/2);e1=n*(n+1); mm=4*m-1;t=(pi/(4*n+2))*(3: 4: mm);nn=(1-(1-1/n)/(8*n*n)); xo=nn*cos(t); forj=1: iter pkm1=1;pk=xo; fork=2: n t1=xo.*pk;pkp1=t1-pkm1-(t1-pkm1)/k+t1; pkm1=pk;pk=pkp1; end den=1.-xo.*xo;d1=n*(pkm1-xo.*pk);dpn=d1./den; d2pn=(2.*xo.*dpn-e1.*pk)./den; d3pn=(4*xo.*d2pn+(2-e1).*dpn)./den; d4pn=(6*xo.*d3pn+(6-e1).*d2pn)./den; u=pk./dpn;v=d2pn./dpn; h=-u.*(1+(.5*u).*(v+u.*(v.*v-u.*d3pn./(3*dpn)))); p=pk+h.*(dpn+(.5*h).*(d2pn+(h/3).*(d3pn+.25*h.*d4pn))); dp=dpn+h.*(d2pn+(.5*h).*(d3pn+h.*d4pn/3)); h=h-p./dp;xo=xo+h; end bp=-xo-h; fx=d1-h.*e1.*(pk+(h/2).*(dpn+(h/3).*(... d2pn+(h/4).*(d3pn+(.2*h).*d4pn)))); wf=2*(1-bp.^2)./(fx.*fx); if(m+m)>n,bp(m)=0;end if~((m+m)==n),m=m-1;end jj=1: m;n1j=(n+1-jj);bp(n1j)=-bp(jj);wf(n1j)=wf(jj); %end 在matlab命令窗口中输入: symsx; f=exp(-x); a=0;b=1;n=3; I=question4(f,a,b,n) 运行结果: I=0.632120255664068 五、插值与逼近 1.给定 上的函数 ,请做如下的插值逼近: (1)构造等距节点分别取 , , 的Lagrange插值多项式; (2)取Chebyshev多项式 的零点: , 作插值节点构造 的插值多项式 和上述的插值多项式均要求画出曲线图形(用不同的线型或颜色表示不同的曲线)。 (1)Lagrange插值matlab程序如下: functionf=L1(n) symsts; q=1/(1+s^2); fori=1: n x(i)=-1+2*(i-1)/(n-1); y(i)=subs(q,s,x(i)); end f=0.0; for(i=1: n) l=y(i); for(j=1: i-1) l=l*(t-x(j))/(x(i)-x(j)); end; for(j=i+1: n) l=l*(t-x(j))/(x(i)-x(j));%计算拉格朗日基函数 end; f=f+l;%计算拉格朗日插值函数 f=collect(f);%化简 end X=-1: 0.1: 1; Y=subs(q,s,X); F=subs(f,t,X); plot(X,Y,'g',X,F,'r'); f=collect(f);%展开 digits(5); f=vpa(f); 1、在matlab命令窗口中输入: f=L1(5) 运行结果 2、在matlab命令窗口中输入: f=L1(8) 运行结果 3、在matlab命令窗口中输入: f=L1(10) 运行结果: 结果分析: 由图像可以看出n=5时拟合曲线与原始曲线的误差比较大,而n=8,10时拟合曲线与原始曲线几乎重合在一起。 可见节点个越多,拉格朗日插值法所得的多项式曲线与原始曲线越接近,误差就越小。 (2)取Chebyshev多项式 的零点为插值点,此时Lagrange插值matlab程序主体基本不变,变的只是插值点,其程序如下: functionf=L2(n) symsts; q=1/(1+s^2); fori=1: n x(i)=cos((2*i-1)*pi/2/n); y(i)=subs(q,s,x(i)); end f=0.0; for(i=1: n) l=y(i); for(j=1: i-1) l=l*(t-x(j))/(x(i)-x(j)); end; for(j=i+1: n) l=l*(t-x(j))/(x(i)-x(j));%计算拉格朗日基函数 end; f=f+l;%计算拉格朗日插值函数 end X=-1: 0.1: 1;%绘图 Y=subs(q,s,X); F=subs(f,t,X); plot(X,Y,'k',X,F,'rp'); f=collect(f);%展开 digits(5); f=vpa(f); 在matlab命令窗口中输入: f=L2(10) 运行结果: 结果分析: 有图像比较可以看出,N=10时,chebyshev零为拉格朗日插值时,其拟合曲线与原始曲线几乎重合,精度较高。 2.已知函数值 3 4 6 6 0 2 和边界条件: . 求三次样条插值函数 并画出其图形. 三次样条插值matlab程序如下: functionThrSample1(x,y,y_1,y_N) symst; f=0.0; if(length(x)==length(y)) n=length(x); else disp('x和y的维数不相等! '); return; end%维数检查 A=diag(2*ones(1,n));%求解m的系数矩阵 u=zeros(n-2,1); lamda=zeros(n-1,1); c=zeros(n,1); fori=2: n-1 u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1)); lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1)); c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+... 3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i)); A(i,i+1)=u(i-1); A(i,i-1)=lamda(i);%形成系数矩阵及向量c end c (1)=2*y_1; c(n)=2*y_N; m=followup(A,c);%用追赶法求解方程组 j=1; whilej<=n-1 h=x(j+1)-x(j); f=y(j)*(2*(t-x(j))+h)*(t-x(j+1))^2/h^3+... y(j+1)*(2*(x(j+1)-t)+h)*(t-x(j))^2/h^3+... m(j)*(t-x(j))*(x(j+1)-t)^2/h^2-... m(j+1)*(x(j+1)-t)*(t-x(j))^2/h^2; f=collect(f); fprintf('当x属于[%d,%d],f=\n',x(j),x(j+1)); disp(f); X=x(j): 0.1: x(j+1);%绘图 F=subs(f,t,X); plot(X,F,'g'); holdon; j=j+1; end functionx=followup(A,b) n=rank(A); for(i=1: n) if(A(i,i)==0) disp('Error: 对角有元素为0! '); return; end end; d=ones(n,1); a=ones(n-1,1); c=ones(n-1); for(i=1: n-1) a(i,1)=A(i+1,i); c(i,1)=A(i,i+1); d(i,1)=A(i,i); end d(n,1)=A(n,n); for(i=2: n) d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1); b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1); end x(n,1)=b(n,1)/d(n,1); for(i=(n-1): -1: 1) x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1); end 在matlab命令窗口中输入: x=[346]; y=[602]; y_1=1; y_N=-1; ThrSample1(x,y,y_1,y_N) 运行结果: 结果分析: 由图像可以看出三次样条插值曲线的优越性,即光滑性好,收敛得到保证,局部性好。 3.观察物体的直线运动,得到如下数据 时刻t 0 0.9 1.9 3.0 3.9 5.0 位移s 0 10 30 51 80 111 求运动学方程 。 matlab程序如下: t=[00.91.93.03.95.0]; s=[010305180111]; n=length(t); a2=0; fori=1: n a2=t(i)^2+a2; end a3=0; fori=1: n a3=t(i)^3+a3; end a4=0; fori=1: n a4=t(i)^4+a4; end b1=0; fori=1: n b1=s(i).*t(i)+b1; end b2=0; fori=1: n b2=s(i).*(t(i).^2)+b2; end A=[a2a3;a3a4] b=[b1;b2] c=inv(A)*b x=0: 0.01: 5.0; y=c (1).*x+c(
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 数值 分析