矩阵与数值分析上机作业满分Word格式文档下载.docx
- 文档编号:22212700
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:36
- 大小:189.64KB
矩阵与数值分析上机作业满分Word格式文档下载.docx
《矩阵与数值分析上机作业满分Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《矩阵与数值分析上机作业满分Word格式文档下载.docx(36页珍藏版)》请在冰豆网上搜索。
x的1范数是5.19
x的2范数是1.28
y的1范数是5050.00
y的2范数是581.68
y的无穷范数是100.00
inputN=1000
x的1范数是7.49
y的1范数是500500.00
y的2范数是18271.11
y的无穷范数是1000.00
第二题
考虑
其中定义f(0)=1,此时f(x)是连续函数。
用此公式计算当
时的函数值,画出图像。
另一方面,考虑下面算法:
d=1+x
ifd=1then
y=1
else
y=lnd/(d-1)
endif
用此算法计算
比较一下发生了什么?
functionproblem2()
N=1000;
n=2*10^(-15)/N;
%公式计算
f=zeros(1,N+1);
t=1;
fori=-10^(-15):
n:
10^(-15)
if(i~=0)
f(t)=log(1+i)/i;
else
f(t)=1;
end
t=t+1;
%算法计算
a=zeros(1,N+1);
d=1+i;
if(d~=1)
a(t)=log(d)/(d-1);
a(t)=1;
%画图,左边是公式算的,右边是算法算的
subplot(1,2,1);
plot(-10^(-15):
10^(-15),f);
holdon
subplot(1,2,2);
10^(-15),a);
第三题
首先编制一个利用秦九韶算法计算一个多项式在给定点的函数值的通用程序,你的程序包括输入多项式的系数以和给定点,输出函数值。
利用你编制的程序计算
在x=2邻域附近的值。
画出
上的图像。
functionproblem3()
%秦九韶算法
n=9;
c=[1,-18,144,-672,2016,-4032,5376,-4608,2304,-512];
N=100;
m=(2.05-1.95)/N;
y=zeros(1,N+1);
ydex=0;
forx=1.95:
m:
2.05
ydex=ydex+1;
dex=1;
temp=c(dex);
whiledex<
=n
dex=dex+1;
temp=temp*x+c(dex);
y(ydex)=temp;
plot(1.95:
2.05,y);
第四题
编制计算给定矩阵A的LU分解和PLU分解的通用程序,然后用你编制的程序完成下面两个计算任务:
(1)考虑
自己取定
,并计算b=Ax。
然后用你编制的不选主元和列主元的Gauss消去法求解该方程组,记你计算出的解为
对n从5到30估计计算解的精度。
(2)对n从5到30计算其逆矩阵。
%LU分解
functionproblem4()
%生成矩阵A
forn=5:
30
A=-ones(n);
A(:
n)=ones(n,1);
foria=1:
n
A(ia,ia)=1;
forja=(ia+1):
(n-1)
A(ia,ja)=0;
x=rand(n,1);
b=A*x;
%LU分解
[n,n]=size(A);
L=eye(n);
M=eye(n);
n-1
ifA(i,i)==0
fprintf('
Error:
A(%d,%d)=0!
\n'
i,i);
return;
forj=i+1:
L(j,i)=A(j,i)/A(i,i);
fork=i+1:
forl=i:
A(k,l)=A(k,l)-L(k,i)*A(i,l);
U=A;
M=inv(U)*inv(L)*M
y1=inv(L)*b;
x1=inv(U)*y1;
error=norm(x-x1)
return;
%PLU分解
%生成矩阵A
n=5;
Q=eye(n);
R=eye(n);
S=eye(n);
L=zeros(n,n,n-1);
foril=1:
L(:
:
il)=eye(n);
P=zeros(n,n,n-1);
P(:
m=i;
a=abs(A(i,i));
ifabs(A(k,i))>
a
a=abs(A(k,i));
m=k;
A([i,m],:
)=A([m,i],:
P([i,m],:
i)=P([m,i],:
i);
L(j,i,i)=-A(j,i)/A(i,i);
A(j,l)=A(j,l)+L(j,i,i)*A(i,l);
U
forx=1:
Q=L(:
x);
fory=x+1:
Q=P(:
y)*Q*P(:
y);
S=S*inv(Q);
%L矩阵
S
forz=1:
R=P(:
z)*R;
%P矩阵
R
b1=R*b;
b2=inv(S)*b1;
x2=inv(U)*b2;
M=inv(U)*inv(S)*R
error=norm(x-x2)
由于篇幅有限仅列出n=5时的精度和逆矩阵
LU分解
M=
0.5000-0.2500-0.1250-0.0625-0.0625
00.5000-0.2500-0.1250-0.1250
000.5000-0.2500-0.2500
0000.5000-0.5000
0.50000.25000.12500.06250.0625
error=
8.8818e-16
PLU分解
7.8516
第五题
编制计算对称正定阵的Cholesky分解的通用程序,并用你编制的程序计算Ax=b,其中
b可以由你自己取定,对n从10到20验证程序的可靠性。
functionproblem5()
%cholesky分解
N=10;
A=zeros(N);
forja=1:
A(ia,ja)=1/(ia+ja-1);
L=chol(A);
L
由于篇幅有限仅列出n=10时的L矩阵
L=
1.00000.50000.33330.25000.20000.16670.14290.12500.11110.1000
00.28870.28870.25980.23090.20620.18560.16840.15400.1417
000.07450.11180.12780.13310.13310.13040.12650.1220
0000.01890.03780.05250.06300.07020.07480.0777
00000.00480.01190.01950.02650.03260.0378
000000.00120.00360.00680.01030.0139
0000000.00030.00110.00220.0038
00000000.00010.00030.0007
000000000.00000.0001
0000000000.0000
第六题
(1)编制程序House(x),其作用是对输入的向量x,输出单位向量u使得
(2)编制Householder变换阵
乘以
的程序HA,注意,你的程序并不显式的计算出H。
(3)考虑矩阵
用你编制的程序计算H使得HA的第一列为
的形式,并将HA的结果显示。
functionproblem6()
x=[1;
-1;
-2;
-10^(1/2);
0];
ifx
(1)>
sg=1;
else
sg=-1;
Nx=length(x);
e1=zeros(Nx,1);
e1
(1)=1;
sum=0;
Nx
sum=sum+x(i)*x(i);
x2=sqrt(sum);
clearsum;
w=x-sg*x2*e1;
u=w/sqrt(w'
*w);
u
A=[1,2,3,4;
-1,3,2^(1/2),3^(1/2);
-2,-2,exp
(1),pi;
-10^(1/2),2,-3,7;
0,2,7,5/2];
H=eye(Nx)-2*u*u'
;
H*A
u=
-0.6124
-0.2041
-0.4082
-0.6455
0
ans=
4.0000-0.83111.4090-6.5378
0.00002.05630.8839-1.7805
0.0000-3.88741.6576-3.8836
0.0000-0.9843-4.6770-4.1078
02.00007.00002.5000
第七题
用Jacobi和Gauss-Seidel迭代求解下面的方程组,输出迭代每一步的误差
:
functionproblem7()
inter=10;
%迭代次数
rx=[69*7/61-8,69*3/(61*2)-7/2,69/61]'
%根据题意生成A,b
A=[5-1-3;
-124;
-3415];
b=[-2110]'
N=length(b);
x=ones(N,1);
%Jacobi迭代法
L=zeros(N,N);
D=L;
U=L;
D(ia,ia)=A(ia,ia);
ifia<
foriu=(ia+1):
U(ia,iu)=A(ia,iu);
ifia>
1
foril=1:
(ia-1)
L(ia,il)=-A(ia,il);
forintt=1:
inter
x=inv(D)*(L+U)*x+inv(D)*b;
wuca1=0;
%求误差
forix=1:
wuca1=wuca1+(x(ix)-rx(ix))^2;
wuca1
%G-S迭代法
x=inv(D-L)*U*x+inv(D-L)*b;
wuca2=0;
wuca2=wuca2+(x(ix)-rx(ix))^2;
wuca2
结果:
wuca1=
24.6036
12.3259
1.9276
5.4520
15.9787
10.0748
3.6731
6.2770
11.9670
8.9591
wuca2=
5.5096
9.1425
6.9372
8.0757
7.4675
7.7807
7.6170
7.7018
7.6577
7.6806
第八题
取不同的初值用Newton迭代以和弦截法求方程
的实根,列表或者画图说明收敛速度。
functionproblem8()
%inNew记录Newton迭代法的收敛速度
%inS记录弦截法的收敛次数
chuzhi=4;
%初值
f=[];
%Newton迭代
x=chuzhi;
forinNew=1:
xx=x-(x^3+2*x^2+10*x-100)/(3*x^2+4*x+10);
f(inNew)=x;
ifabs(xx-x)<
10^(-5)
break;
x=xx;
plot(1:
inNew,f);
clearf;
inNew
xx
%弦截法
xk=chuzhi;
xkr=chuzhi-1;
forinS=1:
xka=xk-(xk^3+2*xk^2+10*xk-100)*(xk-xkr)/((xk^3+2*xk^2+10*xk-100)-(xkr^3+2*xkr^2+10*xkr-100));
f(inS)=xk;
ifabs(xka-xk)<
10^(-5)%精度10^(-5)
xkr=xk;
xk=xka;
inS,f);
inS
xka
inNew=
7
xx=
3.4606
inS=
10
xka=
第九题
用二分法求方程
在区间
上的所有根。
functionproblem9()
%二分法求解方程
%查看根的个数
N=50;
yk=zeros(1,N+1);
h=4*pi/N;
e=[];
c=[];
i=1;
forti=0:
h:
4*pi
e(i)=-2/(exp(ti));
c(i)=cos(ti);
i=i+1;
plot(0:
4*pi,e);
4*pi,c);
%根的取值范围
min=[14710];
max=[25811];
accr=zeros(1,4);
acc=0.0001;
%精度
forr=1:
4
while(max(r)-min(r)>
acc)
ymax=exp(max(r))*cos(max(r))+2;
temp=(max(r)+min(r))/2;
ytemp=exp(temp)*cos(temp)+2;
ifytemp*ymax==0
ifytemp*ymax>
max(r)=temp;
min(r)=temp;
accr(r)=temp;
accr
accr=
1.88074.69407.854810.9955
第十题
考虑函数
用等距节点作
的Newton插值,画出插值多项式以和
的图像,观察收敛性。
functionf=Newton(x,y)
symst;
n=length(x);
c(1:
n)=0;
f=y
(1);
p=1;
y1(j)=(y(j)-y(i))/(x(j)-x(i));
c(i)=y1(i+1);
p=p*(t-x(i));
f=f+c(i)*p;
y=y1;
f=simplify(f);
所以考虑函数f(x)=sin(πx),x∈[0,1],取二次、三次、四次插值函数。
二次:
>
x=linspace(0,1,3)
y=sin(pi*x);
f2=Newton(x,y)
f2=-4*t*(t-1)
三次:
x=linspace(0,1,4);
f3=Newton(x,y)
f3=-(9*3^(1/2)*t*(t-1))/4
四次:
x=linspace(0,1,5);
f4=Newton(x,y)
f4=
(t*(32968342914557536*t^3-65936685829115088*t^2+5181631624232978*t+481984*2^(1/2)+23187))/90992
第十一题
对函数
,取不同的节点数n,用等距节点lagrange插值,观察Runge现象。
functionp=Lagrange(n,t)
x=linspace(-5,5,n);
y=1./(1+x.^2);
m=length(t);
m
s=0;
fork=1:
l(k)=1;
forj=1:
ifj~=k
l(k)=l(k)*(t(i)-x(j))/(x(k)-x(j));
s=l(k)*y(k)+s;
p(i)=s;
end
t=-5:
0.1:
5;
p=1./(1+x.^2);
plot(x,p)
holdon
p5=Lagrange(5,x);
plot(x,p5,'
r'
)
p10=Lagrange(10,x);
plot(x,p10,'
g'
p15=Lagrange(15,x);
plot(x,p15,'
m'
p20=Lagrange(20,x);
plot(x,p20,'
y'
第十二题
令
,考虑积分
区间分为50,100,200,500,1000等,分别用复合梯形以和复合Simpson积分公式计算积分值,将数值积分的结果与精确值比较,列表说明误差的收敛性。
functionproblen12()
请输入N='
%改变区间份数
h=2*pi/N;
%复化梯形公式
x=0:
2*pi;
forin=2:
sum=sum+exp(3*x(in))*cos(pi*x(in));
Tn=2*pi*(1+2*sum+exp(6*pi)*cos(pi*2*pi))/(2*N)
%Simapson
sum1=0;
sum2=sum;
forsi=1:
xx=(h/2)+(si-1)*h;
sum1=sum1+exp(3*xx)*cos(pi*xx);
Sn=2*pi*(1+4*sum1+2*sum2+exp(6*pi)*cos(pi*2*pi))/(6*N)
请输入N=50
Tn=
3.5125e+07
Sn=
3.5231e+07
请输入N=100
3.5205e+07
3.5232e+07
请输入N=200
3.5226e+07
请输入N=500
请输入N=1000
第十三题
分别用2点,3点以和5点的Gauss型积分公式计算如下定积分:
(1)
(2)
functiony=Chebyshev(f,n)
%使用Chebyshev多项式进行Gauss型数值积分
%输入:
f为表达式,n为Gauss点数
%输出:
y为积分值
x=zeros(1,n);
A=zeros(1,n);
x(i)=cos((2*n-1)*pi/(2*n));
A(i)=pi/n;
y=sum(A(1:
n).*subs(f,x(1:
n)));
y=double(y);
functiony=Lengder(f,a,b,x,w)
%使用Lengder
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 矩阵 数值 分析 上机 作业 满分