大连理工大学矩阵大作业.docx
- 文档编号:30500586
- 上传时间:2023-08-16
- 格式:DOCX
- 页数:36
- 大小:541.17KB
大连理工大学矩阵大作业.docx
《大连理工大学矩阵大作业.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵大作业.docx(36页珍藏版)》请在冰豆网上搜索。
大连理工大学矩阵大作业
(1)从大到小的顺序的计算程序:
functiony=snd(n)
formatlong
y=0;
ifn<2
disp('请输入大于1的数!
')
else
s=0;
i=2;
whilei<=n
s=single(s+(1/(i^2-1)));
i=i+1;
end
y=s;
end
(2)从小到大的顺序的计算程序:
functiony=snx(n)
formatlong
y=0;
ifn<2
disp('请输入大于1的数!
')
else
s=0;
i=n;
while1
s=single(s+(1/(i^2-1)));
i=i-1;
ifi==1
break
end
end
y=s;
end
(3)按两种顺序分别计算
并指出有效位数(编制程序时用单精度)
1
的计算结果:
2
的计算结果:
3
的计算结果:
计算时的有效位数为七位数。
1 秦九昭算法计算程序:
functiony=qjz(a,x)
j=3;
i=size(a,2);
switchi
case1
y=a
(1);
case2
y=a
(1)*x+a
(2);
otherwise
p=a
(1)*x+a
(2);
whilej<=i
p=p*x+a(j);
j=j+1;
end
y=p;
end
2 计算
在点23的值。
计算结果如下:
当
时
。
1 Gauss法计算程序和结果:
程序:
A(1,:
)=[31,-13,0,0,0,-10,0,0,0];
A(2,:
)=[-13,35,-9,0,-11,0,0,0,0];
A(3,:
)=[0,-9,31,-10,0,0,0,0,0];
A(4,:
)=[0,0,-10,79,-30,0,0,0,-9];
A(5,:
)=[0,0,0,-30,57,-7,0,-5,0];
A(6,:
)=[0,0,0,0,-7,47,-30,0,0];
A(7,:
)=[0,0,0,0,0,-30,41,0,0];
A(8,:
)=[0,0,0,0,-5,0,0,27,-2];
A(9,:
)=[0,0,0,-9,0,0,0,-2,29];
B=[-15;27;-23;0;-20;12;-7;7;10];
[a,b]=gauss(A,B);
j=size(a,2);
whilej>=1
k=j+1;
s=b(j);
whilek<=9
s=s-x(k)*a(j,k);
k=k+1;
end
x(j)=s/a(j,j);
j=j-1;
end
disp(x)
function[x,y]=gauss(a,b)
num_i=size(a,1);
j=1;
whilej<=(num_i-1)
i=j+1;
whilei<=num_i
r=a(i,j)/a(j,j);
a(i,:
)=a(i,:
)-r*a(j,:
);
b(i,:
)=b(i,:
)-r*b(j,:
);
i=i+1;
end
j=j+1;
end
x=a;
y=b;
运行的结果为:
。
2 列主元消去法计算程序和结果:
A(1,:
)=[31,-13,0,0,0,-10,0,0,0];
A(2,:
)=[-13,35,-9,0,-11,0,0,0,0];
A(3,:
)=[0,-9,31,-10,0,0,0,0,0];
A(4,:
)=[0,0,-10,79,-30,0,0,0,-9];
A(5,:
)=[0,0,0,-30,57,-7,0,-5,0];
A(6,:
)=[0,0,0,0,-7,47,-30,0,0];
A(7,:
)=[0,0,0,0,0,-30,41,0,0];
A(8,:
)=[0,0,0,0,-5,0,0,27,-2];
A(9,:
)=[0,0,0,-9,0,0,0,-2,29];
B=[-15;27;-23;0;-20;12;-7;7;10];
[a,b]=lzy(A,B);
j=size(a,2);
whilej>=1
k=j+1;
s=b(j);
whilek<=9
s=s-x(k)*a(j,k);
k=k+1;
end
x(j)=s/a(j,j);
j=j-1;
end
disp(x)
function[a1,b1]=lzy(a,b)
[num_i,num_j]=size(a);
ab=zeros(num_i,num_j+1);
fork=1:
num_j
ab(:
k)=a(:
k);
end
ab(:
num_j+1)=b(:
1);
j=1;
whilej [max,max_i]=searmax(j,j,ab); i_nu=ab(j,: ); ab(j,: )=ab(max_i,: ); ab(max_i,: )=i_nu; m=j+1; whilem<=num_i forn=j: num_j+1 ab(m,n)=ab(m,n)-(ab(m,j)/max)*ab(j,n); end m=m+1; end j=j+1; end a1=zeros(num_i,num_j); fork=1: num_i a1(: k)=ab(: k); end b1=ab(: num_i+1); function[b,c]=searmax(i,j,a) num_i=size(a,1); k=i; m=abs(a(k,j)); c=k; whilek k=k+1; ifm>=abs(a(k,j)) continue else m=abs(a(k,j)); c=k; end end b=a(c,j); 运行的结果为: (1)LU分解的计算程序及结果: function[l,u]=lufz(a) [num_i,num_j]=size(a); l=eye(num_i); u=eye(num_i); forj=1: num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2; whilei<=num_i j=i; whilej s=0; fork=1: i-1 s=s+l(i,k)*u(k,j); end u(i,j)=a(i,j)-s; s=0; fork=1: i-1 s=s+l(j+1,k)*u(k,i); end l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; fork=1: i-1 s=s+l(i,k)*u(k,num_i); end u(i,num_i)=a(i,num_i)-s; i=i+1; end 输入要求解的矩阵后运行的结果为: (2)带列主元的LU分解计算程序及结果 由于Matlab中自带带列主元的LU分解函数,故计算程序如下: A(1,: )=[31,-13,0,0,0,-10,0,0,0]; A(2,: )=[-13,35,-9,0,-11,0,0,0,0]; A(3,: )=[0,-9,31,-10,0,0,0,0,0]; A(4,: )=[0,0,-10,79,-30,0,0,0,-9]; A(5,: )=[0,0,0,-30,57,-7,0,-5,0]; A(6,: )=[0,0,0,0,-7,47,-30,0,0]; A(7,: )=[0,0,0,0,0,-30,41,0,0]; A(8,: )=[0,0,0,0,-5,0,0,27,-2]; A(9,: )=[0,0,0,-9,0,0,0,-2,29]; [L,U,P]=lu(A); 运行结果如下: 为单位阵。 (3)求A的逆矩阵的计算程序及结果: A(1,: )=[31,-13,0,0,0,-10,0,0,0]; A(2,: )=[-13,35,-9,0,-11,0,0,0,0]; A(3,: )=[0,-9,31,-10,0,0,0,0,0]; A(4,: )=[0,0,-10,79,-30,0,0,0,-9]; A(5,: )=[0,0,0,-30,57,-7,0,-5,0]; A(6,: )=[0,0,0,0,-7,47,-30,0,0]; A(7,: )=[0,0,0,0,0,-30,41,0,0]; A(8,: )=[0,0,0,0,-5,0,0,27,-2]; A(9,: )=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); A_n=eye(num_i); E=eye(num_i); [L,U]=lufz(A); fori=1: num_i y=hd(1,L,E(: i)); A_n(: i)=hd(0,U,y); end disp(A_n) function[l,u]=lufz(a) [num_i,num_j]=size(a); l=eye(num_i); u=eye(num_i); forj=1: num_j u(1,j)=a(1,j); l(j,1)=a(j,1)/u(1,1); end i=2; whilei<=num_i j=i; whilej s=0; fork=1: i-1 s=s+l(i,k)*u(k,j); end u(i,j)=a(i,j)-s; s=0; fork=1: i-1 s=s+l(j+1,k)*u(k,i); end l(j+1,i)=(a(j+1,i)-s)/u(i,i); j=j+1; end s=0; fork=1: i-1 s=s+l(i,k)*u(k,num_i); end u(i,num_i)=a(i,num_i)-s; i=i+1; end functionx=hd(f,a,b) [num_i,num_j]=size(a); x=zeros(num_i,1); switchf case0 i=num_i-1; x(num_i)=b(num_i)/a(num_i,num_j); whilei>=1 s=0; fork=i+1: num_i s=s+a(i,k)*x(k); end x(i)=(b(i)-s)/a(i,i); i=i-1; end case1 x (1)=b (1)/a(1,1); j=2; whilej<=num_i s=0; fork=1: j-1 s=s+a(j,k)*x(k); end x(j)=(b(j)-s)/a(j,j); j=j+1; end otherwise disp('error! 请输入正确的参数! ') end 运行结果: (4)求A的行列式的计算程序及结果: A(1,: )=[31,-13,0,0,0,-10,0,0,0]; A(2,: )=[-13,35,-9,0,-11,0,0,0,0]; A(3,: )=[0,-9,31,-10,0,0,0,0,0]; A(4,: )=[0,0,-10,79,-30,0,0,0,-9]; A(5,: )=[0,0,0,-30,57,-7,0,-5,0]; A(6,: )=[0,0,0,0,-7,47,-30,0,0]; A(7,: )=[0,0,0,0,0,-30,41,0,0]; A(8,: )=[0,0,0,0,-5,0,0,27,-2]; A(9,: )=[0,0,0,-9,0,0,0,-2,29]; [num_i,num_j]=size(A); [L,U]=lufz(A); s=1; fori=1: num_i s=s*U(i,i); end disp(['矩阵的行列式值为: ',num2str(s)]) 运行程序后结果与调用matlab中det()函数结果如下: (1)求cholesky分解程序及结果: functionl=chole(a) [num_i,num_j]=size(a); l=zeros(num_i); l(1,1)=sqrt(a(1,1)); fork=2: num_i l(k,1)=a(k,1)/l(1,1); end j=2; whilej<=num_j s1=0; fork=1: j-1 s1=s1+l(j,k)^2; end l(j,j)=sqrt(a(j,j)-s1); i=j+1; whilei<=num_i s2=0; fork=1: j-1 s2=s2+l(i,k)*l(j,k); end l(i,j)=(a(i,j)-s2)/l(j,j); i=i+1; end j=j+1; end 运行程序后的结果为: (2)求解方程组程序及结果: a=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11]; b=[13;-9;6;0]; L=chole(a); y=hd(1,L,b); x=hd(0,L',y); disp(x) 运行后的结果为: 程序和结果如下: A=[4,2,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6]; B=[6;2;10;5]; [L,U]=sjfj(A);% y=hd(1,L,B); x=hd(0,U,y); disp('方程组的解为: ') disp(x) function[l,u]=sjfj(a) num_i=size(a,1); i=2; whilei<=num_i a(i,i-1)=a(i,i-1)/a(i-1,i-1); a(i,i)=a(i,i)-a(i,i-1)*a(i-1,i); i=i+1; end l=tril(a); u=triu(a); forj=1: num_i; l(j,j)=1; end 运行程序后结果为: 编程求解矩阵A的QR分解: (1)QR分解计算程序: function[Q,R]=hqr(A) [n,n]=size(A); Q=eye(n); fori=1: (n-1) E=eye(n-i+1); e1=E(: 1); a=zeros(1,n-i+1)'; forj=1: (n-i+1) a(j)=A(j+i-1,i); end av=sqrt(a'*a); w=a-av*e1; h=E-(2/(w'*w))*(w*w'); q=eye(n); fork=i: n forj=i: n q(k,j)=h(k-i+1,j-i+1); end end A=q*A; Q=q*Q; end R=A; Q=Q'; (2)输入矩阵A后的计算结果: (1)Jacobi迭代法求解程序与结果: a=0; b=0; c=0; while1 a0=0.25*(7+b-c); b0=(1/8)*(-21-4*a-c); c0=(1/5)*(15+2*a-b); A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); iff<=0.001 x=[a0,b0,c0]; break else a=a0; b=b0; c=c0; end end disp('方程组的解为x') disp(x) 运行后的结果为: (2)Gauss-Seidel迭代法求解的程序与结果: a=0; b=0; c=0; while1 a0=a; b0=b; c0=c; a=(1/4)*(7+b-c); b=(1/8)*(-21-4*a-c); c=(1/5)*(15+2*a-b); A=[abs(a-a0),abs(b-b0),abs(c-c0)]; f=max(A); iff<=0.001 x=[a,b,c]; break end end disp('方程组的解为: ') disp(x) 运行后的结果为: (1)计算 在 上的根的程序: symsx f=2*x^3-5*x-1; df=diff(f,x); g=x-(f/df); x0=1; while1 x1=subs(g,x0); dx=abs(x1-x0); ifdx<=0.001 disp(['Newton法的解为: x=',num2str(x1)]) break else x0=x1; end end x0=1; x1=1.1; while1 x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); ifdx<=0.001 disp(['割线法的解为: x=',num2str(x2)]) break else x0=x1; x1=x2; end end 运行后结果为: Newton法 ,割线法 。 (2)计算 在 上的根的程序: symsx f=exp(x)*sin(x); df=diff(f,x); g=x-(f/df); x0=-2; while1 x1=subs(g,x0); dx=abs(x1-x0); ifdx<=0.001 disp(['Newton法的解为: x=',num2str(x1)]) break else x0=x1; end end x0=-2; x1=-2.1; while1 x2=x1-(subs(f,x1)/(subs(f,x1)-subs(f,x0)))*(x1-x0); dx=abs(x2-x1); ifdx<=0.001 disp(['割线法的解为: x=',num2str(x2)]) break else x0=x1; x1=x2; end end 运行后结果为: Newton法 ,割线法 。 symsx f=x*cos(x)-2; x0=-4; x2=-2; while1 x1=0.5*(x0+x2); f0=subs(f,x0); f1=subs(f,x1); c=f0*f1; ifc<0 x2=x1; else x0=x1; end l=abs(x2-x0); ifl<=0.001 disp(['二分法的解为x=',num2str(x0)]) break end end 运行后的结果为 。 (1)Lagrange插值程序: functionyv=lagran(xd,yd,xv) symsx num=length(xd);%计算节点的个数 f=0; fork=1: num%创建lagrange插值函数 index=setdiff(1: num,k); f=f+yd(k)*prod((x-xd(index))./(xd(k)-xd(index))); end yv=subs(f,xv);%计算xv点插值函数值 (2)绘制函数图程序: symsx f=1./(1+x.^2); xd2=-5: 2: 5; yd2=subs(f,xd2); x0=-5: 0.1: 5; y=subs(f,x0); plot(x0,y) holdon y2=lagran(xd2,yd2,x0); plot(x0,y2) (3)运行程序后得到的步长分别为2,1, 的插值函数图与原图形的比较图如下: (1)三次样条插值程序: functionyv=csi(xd,yd,xv,h,mf,ml) num_xd=length(xd); n=num_xd-1; a_m=1: num_xd-2; fori=1: num_xd-2 a_m(i)=2; end a_np=1: num_xd-3; fori=1: num_xd-3 a_np(i)=0.5; end A=diag(a_m)+diag(a_np,-1)+diag(a_np,1); g=1: n-1; fork=1: n-1 g(k)=3*((0.5/h)*(yd(k+2)-yd(k+1))+(0.5/h)*(yd(k+1)-yd(k))); end b=1: n-1; b (1)=g (1)-0.5*mf; b(n-1)=g(n-1)-0.5*ml; fork=2: n-2 b(k)=g(k); end b=b.'; M=A\b; m=1: n+1; m (1)=mf; m(n+1)=ml; fori=2: n m(i)=M(i-1); end num_xv=length(xv); yv=1: num_xv; fori=1: num_xv forj=1: num_xd ifxv(i) k1=j-1; k2=k1+1; break end end yv1=(h+2*(xv(i)-xd(k1)))*(xv(i)-xd(k2))^2*(yd(k1)/h^3); yv2=(h-2*(xv(i)-xd(k2)))*(xv(i)-xd(k1))^2*(yd(k2)/h^3); yv3=(xv(i)-xd(k1))*(xv(i)-xd(k2))^2*(m(k1)/h^2); yv4=(xv(i)-xd(k2))*(xv(i)-xd(k1))^2*(m(k2)/h^2); yv(i)=yv1+yv2+yv3+yv4; end (2)绘图程序(改变程序中的h值即可改变步长): clearall clc h=0.5; xd=-5: h: 5; yd=1./(1+xd.^2); xv=-5: 0.1: 5; y=1./(1+xv.^2); mf=0.014793; ml=-0.014793; yv=csi(xd,yd,xv,h,mf,ml); plot(xv,y) holdon plot(xv,yv) (3)步长为2,1, 时原函数图与插值函数图的比较: (1)复化梯形公式计算积分程序: functiont=trape(f,a,b,n) h=(b-a)/n; index=(a+h): h: (b-h); t1=sum(subs(f,index)); t=((b-a)/(2*n))*(subs(f,a)+2*t1+subs(f,b)); (2)复化Simpson公式计算程序: functions=simpson(f,a,b,n) h=(b-a)/n; index1=(a+0.5*h): h: (b-0.5*h); index2=(a+h): h: (b-h); s1=sum(subs(f,index1)); s2=sum(subs(f
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 大连理工大学 矩阵 作业
![提示](https://static.bdocx.com/images/bang_tan.gif)