矩阵上机题副1.docx
- 文档编号:23220535
- 上传时间:2023-05-15
- 格式:DOCX
- 页数:36
- 大小:782.26KB
矩阵上机题副1.docx
《矩阵上机题副1.docx》由会员分享,可在线阅读,更多相关《矩阵上机题副1.docx(36页珍藏版)》请在冰豆网上搜索。
矩阵上机题副1
矩阵与数值分析上机实习题
王安然(21509201)
工具:
matlab
1.
(1)
(2)
Sn1=single(0);%存储
(1)问的结果
Sn2=single(0);%存储
(2)问的结果
N=10^2;%输入值分别为10^2,10^4,10^6
fori=2:
N
Sn1=Sn1+1/(i^2-1);
end;
fori=N:
-1:
2
Sn2=Sn2+1/(i^2-1);
end;
(3)N分别取102,103,104结果如下
(4)该题说明,改变元素和的相加顺序影响极限值。
计算机中存在大小数相加的舍入误差误差。
方法二的方法要比方法一要好,因为避免了大数吃小数的舍入误差。
2.
x=23;
A=[7;3;-5;11];
m=size(A);
count=(m-1:
-1:
0);%多项式次数向量
f=A'*x.^count';
3.
%initialization
A=[31,-13,0,0,0,-10,0,0,0;
-13,35,-9,0,-11,0,0,0,0;
0,-9,31,-10,0,0,0,0,0;
0,0,-10,79,-30,0,0,0,-9;
0,0,0,-30,57,-7,0,-5,0;
0,0,0,0,-7,47,-30,0,0;
0,0,0,0,0,-30,41,0,0;
0,0,0,0,-5,0,0,27,-2;
0,0,0,-9,0,0,0,-2,29];
b=[-15;27;-23;0;-20;12;-7;7;10];
if(det(A)==0)
disp('NoUniqueSolution');
return;
end;
m=size(A,1);
%GuassianElimination
aug=[A,b];
fori=1:
m
forj=i+1:
m
aug(j,:
)=aug(j,:
)-aug(j,i)/aug(i,i)*aug(i,:
);
end;
end;
disp(aug)
x_GaussElim=solveFun(aug);
fprintf('xforGuassianElimination\n');
disp(x_GaussElim);
%GuassianColumnPCA
aug1=[A,b];
fori=1:
m
max=aug1(i,i);
flag=i;
fork=i+1:
m
if(aug1(k,i)^2>max^2)
max=aug1(k,i);
flag=k;
end;
end;
temp=aug1(flag,:
);
aug1(flag,:
)=aug1(i,:
);
aug1(i,:
)=temp;
forj=i+1:
m
aug1(j,:
)=aug1(j,:
)-aug1(j,i)/aug1(i,i)*aug1(i,:
);
end;
end;
x_GaussColu=solveFun(aug1);
fprintf('xforGaussianColumnPCA\n');
disp(x_GaussColu)
其中解上三角方程组的函数solveFun(aug):
functionx=solveFun(ang)
n=size(ang,1);
x=zeros(n,1);
A=ang(:
1:
n);
Y=ang(:
n+1);
x(n,1)=Y(n,1)/A(n,n);
fori=n-1:
-1:
1
x(i,1)=(Y(i,1)-A(i,i+1:
n)*x(i+1:
n,1))/A(i,i);
end;
disp(x);
4.
%LU
function[L,U]=LU(A)
m=size(A,1);
%LU
A1=A;
L=eye(m);
fori=1:
m
l=eye(m);
forj=i+1:
m
l(j,i)=-A1(j,i)/A1(i,i);
end;
A1=l*A1;
L=L*l^(-1);
end;
U=A1;
fprintf('L:
\n');
disp(L);
fprintf('U:
\n');
disp(U);
end
%A_inverse
functionA_inverse=LU_inverse(A)
[L,U]=LU(A);
n=size(A,1);
L_inv=zeros(n,n);
U_inv=zeros(n,n);
%inv(L)
fori=1:
1:
n
L_inv(i,i)=1/L(i,i);
end
fori=2:
1:
n
forj=1:
1:
i-1
L_inv(i,j)=-L_inv(i,i)*(L(i,j:
(i-1))*L_inv(j:
(i-1),j));
end
end
%inv(U)
fori=1:
1:
n
U_inv(i,i)=1/U(i,i);
end
fori=1:
1:
n-1
k=1;
forj=i+1:
1:
n
U_inv(k,j)=-U_inv(k,k)*(U(k,(k+1):
j)*U_inv((k+1):
j,j));
k=k+1;
end
end
A_inverse=L_inv*U_inv;
fprintf('A_inverse\n');
disp(A_inverse);
end
%A_det
functionA_det=LU_det(A)
m=size(A);
[L,U]=LU(A);
U_det=1;
fori=1:
m
U_det=U_det*U(i,i);
end
A_det=U_det;
fprintf('A_det\n');
disp(A_det);
end
%LUColumnPCA
function[L,U,P]=LU_ColumnPCA(A2)
m=size(A2,1);
P=eye(m);
L=eye(m);
p_cell=cell(1,m-1);
l_cell=cell(1,m-1);
fori=1:
m
l=eye(m);
p=eye(m);
max=A2(i,i);
flag=i;
forj=i+1:
m%getmaxinthecolumn
if(A2(j,i)^2>max^2)%compareabsolutevalue
max=A2(j,i);
flag=j;
end;
end;
temp1=p(flag,:
);
p(flag,:
)=p(i,:
);
p(i,:
)=temp1;
A2=p*A2;
forj=i+1:
m
l(j,i)=-A2(j,i)/A2(i,i);
end;
A2=l*A2;
p_cell{i}=p;
l_cell{i}=l;
P=p*P;
end;
U=A2;
fprintf('P=\n');
disp(P);
fork=1:
m-1
L_temp=l_cell{k};
forn=k+1:
m-1
L_temp=p_cell{n}*L_temp*p_cell{n}^(-1);
end;
L=L*L_temp^(-1);
end;
fprintf('L=\n');
disp(L);
fprintf('U=\n');
disp(U);
5.
A=[2,1,-1,1;1,5,2,7;-1,2,10,1;1,7,1,11];
b=[13,-9,6,0];
disp(A);
ifA~=(A')
error('要求为对称正定阵');
end;
[m,n]=size(A);
u=zeros(n,n);
fori=1:
n
u(i,i)=A(i,i);
fork=1:
i-1
u(i,i)=u(i,i)-u(k,i)^2;
end
u(i,i)=sqrt(u(i,i));
forj=i+1:
n
u(i,j)=A(i,j);
fork=1:
i-1
u(i,j)=u(i,j)-u(k,i)*u(k,j);
end;
u(i,j)=u(i,j)/u(i,i);
end;
end;
L=u';
disp(L)
y=zeros(n,1);
fori=1:
n%计算y
temp=0;
fork=1:
i-1
temp=temp+L(i,k)*y(k);
end;
y(i)=(b(i)-temp)/L(i,i);
end;
x=zeros(n,1);
fori=n:
-1:
1%计算X
temp=0;
fork=i+1:
n
temp=temp+L(k,i)*x(k);
end;
x(i)=(y(i)-temp)/L(i,i);
end;
disp(x);
6.
A=[4,4,0,0;3,-2,1,0;0,2,5,3;0,0,-1,6];
b=[6;2;10;5];
m=size(A,1);
y=zeros(m,1);
x=zeros(m,1);
[L,U]=LU(A);
y
(1)=b
(1);
fori=2:
m
y(i)=b(i)-L(i,i-1)*y(i-1);
end;
x(m)=y(m)/U(m,m);
forj=m-1:
-1:
1
x(j)=(y(j)-A(j,j+1)*x(j+1))/U(j,j);
end;
fprintf('x=\n');
disp(x);
7.
function[Q,R]=QR(A)
m=size(A);
I=eye(m);
Q=eye(m);
R=A;
fori=1:
m-1
q=zeros(m);
w=R(i:
end,i)-norm(R(i:
end,i))*I(i:
end,i);
h=I(i:
end,i:
end)-2/(w'*w)*(w*w');
q(i:
end,i:
end)=h;
if(i>1)
forj=1:
i-1
q(j,j)=1;
end;
end;
R=q*R;
Q=Q*q';
end;
fprintf('R:
\n');
disp(R);
fprintf('Q:
\n');
disp(Q)
8.
Jabobi方法
function[x,i]=Jacobi(A,b)
m=size(A,1);
x=zeros(m,1);
I=eye(m);
D=I.*A;
fprintf('D:
\n');
disp(D)
L=D-tril(A);
fprintf('L:
\n');
disp(L);
U=D-triu(A);
fprintf('U:
\n');
disp(U);
i=1;
while(true)
x1=D^(-1)*(L+U)*x+D^(-1)*b;
if(norm(x1-x,inf)<=1e-4);
break;
end;
x=x1;
i=i+1;
end;
fprintf('x:
\n');
disp(x);
fprintf('iteratetimes:
\n');
disp(i);
GS方法
function[x,i]=GS(A,b)
m=size(A,1);
x=zeros(m,1);
I=eye(m);
D=I.*A;
disp(D)
L=D-tril(A);
disp(L);
U=D-triu(A);
disp(U);
i=1;
while(true)
x1=(D-L)^(-1)*U*x+(D-L)^(-1)*b;
if(norm(x1-x,inf)<=1e-4);
break;
end;
x=x1;
i=i+1;
end;
fprintf('x:
\n');
disp(x);
fprintf('iteratetimes:
\n');
disp(i);
9.
(1)2*x^3-5*x-1=0;
Newton法
function[x1,i]=newton()
symsx;
f(x)=2*x^3-5*x-1;%f
x0=1.5;
i=0;
while(true)
x=x0;
x1=x0-f(x)/eval(diff(f));
if(abs(x0-x1)<=1e-6);
break;
end;
x0=x1;
i=i+1;
end;
fprintf('x:
\n');
disp(vpa(x1,6));
fprintf('iteratetimes:
\n');
disp(i);
弦割法
function[x2,i]=secant(f,x0,x1)
i=0;
while(true)
x2=x1-f(x1)/(f(x1)-f(x0))*(x1-x0);
if(abs(x2-x1)<=1e-6)
break;
end;
x0=x1;
x1=x2;
i=i+1;
end;
fprintf('x:
\n');
disp(x2);
fprintf('iteratetimes:
\n');
disp(i);
(2)exp
(1)^x*sin(x)=0;
Newton法
f(x)=exp
(1)^x*sin(x);
x0=-3.5;
弦割法
10.二分法求f=x*cos(x)-2=0;
定义函数如下:
functionf=myfun1(x)
f=x*cos(x)-2;
二分法函数:
functionx=bisection(f,a,b)
i=0;
e=b-a;
while(e>=1e-6)
x=(a+b)/2;
if(f(x)*f(a)>0)
%disp(f(x));
a=x;
elseif(f(x)*f(a)<0)
b=x;
else
a=x;b=x;
end;
e=e/2;
i=i+1;
end;
fprintf('x:
\n');
disp(x)
fprintf('iteratetimes:
\n');
disp(i)
11.f=1/(1+x^2)
functionyh=Lagrange_chazhi(x,y,xh)
n=length(x);
m=length(xh);
x=x(:
);
y=y(:
);
xh=xh(:
);
yh=zeros(m,1);
c1=ones(1,n-1);
c2=ones(m,1);
fori=1:
n,
xp=x([1:
i-1i+1:
n]);
yh=yh+y(i)*prod((xh*c1-c2*xp')./(c2*(x(i)*c1-xp')),2);
end
调用函数
xh=[-5:
0.001:
5];
y1=1./(1+xh.^2);
subplot(2,2,1);
plot(xh,y1);
gridon;
title('原图');
%step=2
x1=[-5:
2:
5];
y=1./(1+x1.^2);
yh1=Lagrange_chazhi(x1,y,xh);
subplot(2,2,2)
plot(xh,yh1)
gridon;
title('步长2');
%step=1
x2=[-5:
1:
5];
y=1./(1+x2.^2);
yh2=Lagrange_chazhi(x2,y,xh);
subplot(2,2,3)
plot(xh,yh2)
gridon;
title('步长1');
%step=1/2
x3=[-5:
1/2:
5];
y=1./(1+x3.^2);
yh3=Lagrange_chazhi(x3,y,xh);
subplot(2,2,4)
plot(xh,yh3)
gridon;
title('步长1/2');
12.
a=-5;
b=5;
step=2;
n=(b-a)/step;
f_a=(-2*a)/(1+a^2)^2;
f_b=(-2*b)/(1+b^2)^2;
x=
zeros(1,n+1);
y=zeros(1,n+1);
fori=1:
(n+1)
x(i)=a+(i-1)*step;
y(i)=1/(1+x(i)^2);
end
%给定1阶边界条件
pp=csape(x,y,'complete',[f_a,f_b]);
xi=-5:
0.1:
5;yi=ppval(pp,xi);
plot(x,y,'o',xi,yi);
gridon;
return;
13.f=x^2*sin(x);
梯形法
functionT=TiXing_quad(fun,a,b,n)
x=[a:
(b-a)/n:
b];
x=x(:
);
T=0;
fori=1:
n
T=T+fun(x(i))+fun(x(i+1));
end
T=T*(b-a)/(2*n);
Simpson法
functionT=Simpson_quad(fun,a,b,n)
x=[a:
(b-a)/n:
b];
x=x(:
);
T=0;
fori=1:
n
T=T+fun(x(i))+4*fun((x(i)+x(i+1))/2)+fun(x(i+1));
end
T=T*(b-a)/(6*n);
精确计算:
。
。
。
结论:
通过比较可知,梯形法与Simpson法精度。
。
。
。
14.Guass型求积公式
function[y,Ak,xk]=GaussLegendre(fun,a,b,n)
%y积分结果
%Ak求积系数
%xk求积节点
%计算求积节点
symsx;
tol=1e-8;
%tk为求积节点
p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n));
tk=roots(p);
%计算求积系数
Ak=zeros(n+1,1);
fori=1:
n+1
xkt=tk;
xkt(i)=[];
pn=poly(xkt);
fp=@(x)polyval(pn,x)/polyval(pn,tk(i));
Ak(i)=quadl(fp,-1,1,tol);%求积系数
end
%积分变量代换,将[a,b]变换到[-1,1]
xk=(b-a)/2*tk+(b+a)/2;
%检验积分函数fun有效性
fun=fcnchk(fun,'vectorize');
%计算变量代换之后积分函数的值
fx=fun(xk)*(b-a)/2;
%计算积分值
y=sum(Ak.*fx);
end
15.
定义函数
functionf=myfun2(x,u)
f=2/x*u+(x^2)*(exp
(1)^x);
Euler法
functionEuler(f,a,b,h,u0)%f:
f(x,u)函数,a,b区间下上界,h:
步长,U0=U(1;)
n=(b-a)/h;
u=zeros(n,1);
x=a:
h:
b;
u
(1)=u0;
fori=1:
n
u(i+1)=u(i)+h*f(x(i),u(i));
end
plot(x,u,'r');
fprintf('Euler:
\n');
disp(u);
(1)调用方法(步长为0.1时):
f=@myfun2;
Euler(f,1,2,0.1,0)
(2)调用方法(步长为0.05时):
f=@myfun2;
Euler(f,1,2,0.05,0)
(3)调用方法(步长为0.01时):
f=@myfun2;
Euler(f,1,2,0.01,0)
U向量显示略
改进Euler法
functiongaijin_Euler(f,a,b,h,u0)
n=(b-a)/h;
u=zeros(n,1);
x=a:
h:
b;
u
(1)=u0;
fori=1:
n
u(i+1)=u(i)+h/2*(f(x(i),u(i))+f(x(i+1),u(i)+h*f(x(i),u(i))));
end
plot(x,u,'g');
fprintf('gaijin_Euler:
\n');
disp(u);
(1)调用方法(步长为0.1时):
f=@myfun2;
gaijin_Euler(f,1,2,0.1,0)
(2)调用方法(步长为0.05时):
f=@myfun2;
gaijin_Euler(f,1,2,0.05,0)
(3)调用方法(步长为0.01时):
f=@myfun2;
gaijin_Euler(f,1,2,0.01,0)
U向量显示略
Runge-Kutta法
functionRunge_Kutta_4(f,a,b,h,u0)
n=(b-a)/h;
u=zeros(n,1);
x=a:
h:
b;
u
(1)=u0;
fori=1:
n
k1=f(x(i),u(i));
k2=f(x(i)+1/2*h,u(i)+1/2*h*k1);
k3=f(x(i)+1/2*h,u(i)+1/2*h*k2);
k4=f(x(i)+h,u(i)+h*k3);
u(i+1)=u(i)+h/6*(k1+2*k2+2*k3+k4);
end
plot(x,u,'b');
fprintf('Runge_Kutta_4:
\n');
disp(u);
(1)调用方法(步长为0.1时):
f=@myfun2;
Runge_Kutta_4(f,1,2,0.1,0)
(2)调用方法(步长为0.05时):
f=@myfun2;
Runge_Kutta_4(f,1,2,0.05,0)
(3)调用方法(步长为0.01时):
f=@myfun2;
Runge_Kutta_4(f,1,2,0.01,0)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 上机