大连理工大学矩阵与数值分析上机作业Word格式.docx
《大连理工大学矩阵与数值分析上机作业Word格式.docx》由会员分享,可在线阅读,更多相关《大连理工大学矩阵与数值分析上机作业Word格式.docx(24页珍藏版)》请在冰豆网上搜索。
);
holdon
plot(x,y2);
2.2结果
2.3分析
红色的曲线代表未考虑题中算法时的情况,如果考虑题中的算法则数值大小始终为1,这主要是由于大数加小数的原因。
第3题
3.1程序
A=[1-18144-6722016-40325376-46082304-512];
x=1.95:
0.005:
2.05;
length(x);
y1(i)=f(A,x(i));
y2(i)=(x(i)-2)^9;
figure(3);
plot(x,y1);
holdon;
plot(x,y2,'
F.m文件
functiony=f(A,x)
y=A
(1);
fori=2:
length(A);
y=x*y+A(i);
end;
3.2结果
第4题
4.1程序
A=2*eye(n)-tril(ones(n,n),0);
n
A(i,n)=1;
n=length(A);
U=A;
e=eye(n);
n-1
[max_data,max_index]=max(abs(U(i:
n,i)));
e0=eye(n);
max_index=max_index+i-1;
U=e0*U;
e1=eye(n);
forj=i+1:
e1(j,i)=-U(j,i)/U(i,i);
U=e1*U;
P{i}=e0;
%把变换矩阵存到P中
L{i}=e1;
e=e1*e0*e;
fork=1:
n-2
Ldot{k}=L{k};
fori=k+1:
Ldot{k}=P{i}*Ldot{k}*P{i};
Ldot{n-1}=L{n-1};
LL=eye(n);
PP=eye(n);
PP=P{i}*PP;
LL=Ldot{i}*LL;
b=ones(n,2);
b=e*b;
%解方程
x=zeros(n,1);
x(n)=b(n)/U(n,n);
fori=n-1:
-1:
1
x(i)=(b(i)-U(i,:
)*x)/U(i,i);
X=U^-1*e^-1*eye(n);
%计算逆矩阵
AN=X'
;
result2{n-4,1}=AN;
result1{n-4,1}=x;
fprintf('
%d:
\n'
n)
%d'
AN);
4.2结果
n=5
1.0625
-0.875
-0.75
-0.5
-0.0625
0.0625
1.125
0.125
1.25
0.25
1.5
-0.125
-0.25
n=10
同样的方法可以算出n=20,n=30时的结果,这里就不罗列了。
第5题
5.1程序
10至20'
forj=1:
a(i,j)=1/(i+j-1);
forj=1:
sum=0;
fork=1:
j-1
sum=sum+l(j,k)^2;
l(j,j)=sqrt(a(j,j)-sum);
fori=j+1:
sum=sum+l(i,k)*l(j,k);
l(i,j)=(a(i,j)-sum)/l(j,j);
b=ones(n,1);
y=zeros(n,1);
y(n)=b(n)/l(n,n);
y(i)=(b(i)-l(i,:
)*y)/l(i,i);
l=l'
x(n)=y(n)/l(n,n);
x(i)=(y(i)-l(i,:
)*x)/l(i,i);
%d\t'
x);
5.2结果
n=11
n=12
n=13
n=14
n=15
n=16
n=17
n=18
n=19
n=20
-746517.8688
3111493.423
-11884558.85
478355909.6
497329749.5
519377549.9
445748685
378885969.7
341571897.8
1020943960
994971030.6
69820595.08
-354658479.9
1634110905
-77484610111
-80547115920
-82914903907
-72327488690
-62111481037
-60010758946
-1.76915E+11
-1.68795E+11
-1587444197
9875543407
-54993461592
3.06265E+12
3.18327E+12
3.23484E+12
2.86169E+12
2.48062E+12
2.55037E+12
7.40163E+12
6.93865E+12
152********
-1.17236E+11
7.93546E+11
-5.20373E+13
-5.40791E+13
-5.42806E+13
-4.8713E+13
-4.26882E+13
-4.65382E+13
-1.3123E+14
-1.21156E+14
-76184620048
7.35352E+11
-6.11155E+12
4.7524E+14
4.93812E+14
4.89563E+14
4.46792E+14
3.97243E+14
4.57473E+14
1.22E+15
1.11E+15
2.18036E+11
-2.70378E+12
2.80302E+13
-2.62E+15
-2.72E+15
-2.66E+15
-2.48E+15
-2.25E+15
-6.64E+15
-5.98E+15
-3.70513E+11
6.12295E+12
-8.11E+13
9.24E+15
9.59E+15
9.27E+15
8.90E+15
8.30E+15
1.04E+16
2.17E+16
1.94E+16
3.69292E+11
-8.64269E+12
1.51789E+14
-2.17E+16
-2.25E+16
-2.14E+16
-2.07E+16
-2.60E+16
-4.18E+16
-3.74E+16
-1.99261E+11
7.40507E+12
-1.83339E+14
3.40E+16
3.53E+16
3.34E+16
3.52E+16
3.56E+16
4.28E+16
4.06E+16
3.77E+16
44905979430
-3.52275E+12
1.3792E+14
-3.55E+16
-3.68E+16
-3.47E+16
-3.89E+16
-4.12E+16
-4.38E+16
-4.89E+15
-1.01E+16
7.13565E+11
-5.87483E+13
2.35E+16
2.44E+16
2.37E+16
2.73E+16
2.94E+16
2.40E+16
-1.64E+16
-5.28E+15
1.08203E+13
-8.98E+15
-9.29E+15
-1.03E+16
-9.10E+15
-7.78E+15
-4.50E+15
-1.23E+16
-1.58E+16
1.50E+15
1.55E+15
3.06E+15
-2.94E+15
-6.68E+15
1.72E+15
1.63E+16
2.17168E+12
-7.94676E+14
5.08E+15
7.29E+15
-1.45E+15
3.84E+16
3.85E+16
1.58892E+14
-2.50E+15
-2.19E+15
-8.57E+15
-2.64E+16
-2.01E+16
4.80384E+14
-3.1296E+14
1.48E+16
-8.62E+16
-7.98E+16
2.304E+14
-9.01E+15
1.40E+17
1.19E+17
1.99E+15
-8.07E+16
-6.35E+16
1.70E+16
1.09E+16
7.5453E+14
第6题
6.1程序
A=[1234;
-13sqrt
(2)sqrt(3);
-22exp
(1)pi;
-sqrt(10)2-37;
0275/2;
];
U=f61(A(:
2));
HA=f62(U,A);
f.m文件
functionU=f61(x)
e1=eye(length(x),1);
U=x-sign(x
(1))*sqrt(dot(x,x))*e1;
U=U./sqrt(dot(U,U));
functionHA=f62(U,A)
HA=A-2*U*U'
*A;
6.2结果
-2.264911064
5
4.735840869
7.695867546
2.264911064
4.44E-16
-0.321627306
-1.963816738
0.176607376
1.561054583
0.677680956
-0.985670284
-4.157227246
4.536088303
2.176607376
5.842772754
0.036088303
第7题
7.1程序
max=1000;
x(1,:
)=[123];
max
x(i+1,1)=(-2+x(i,2)+3*x(i,3))/5;
x(i+1,2)=(1+x(i,1)-4*x(i,3))/2;
x(i+1,3)=(10+3*x(i,1)-4*x(i,2))/15;
err(i)=sqrt(dot(x(i+1,:
)-x(i,:
),x(i+1,:
)));
if(err(i)<
10^-6)
break;
figure(7);
plot(err);
clearerr;
x(i+1,2)=(1+x(i+1,1)-4*x(i,3))/2;
x(i+1,3)=(10+3*x(i+1,1)-4*x(i+1,2))/15;
plot(err,'
7.2结果
误差越来越小。
第8题
8.1程序
max=100;
x
(1)=1;
x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(3*x(i)^2+4*x(i)+10);
if(abs(x(i+1)-x(i))<
figure(8)
plot(x);
clearx;
x
(1)=0;
x
(2)=1;
x(i+1)=x(i)-(x(i)^3+2*x(i)^2+10*x(i)-100)/(x(i)^3+2*x(i)^2+10*x(i)-(x(i-1)^3+2*x(i-1)^2+10*x(i-1)))*(x(i)-x(i-1));
plot(x,'
8.2结果
8.3分析
由计算结果可知,弦截法的收敛速度比牛顿法的收敛速度快。
第9题
9.1程序
f(0,4*pi);
function[]=f(l,r)
if(r-l<
10^-6)
fprintf('
%g,'
(r-l)/2+l);
return
iff9(l)*f9(r)<
0
iff9((l+r)/2+l)*f9(l)<
f(l,(l+r)/2+l);
iff9((l+r)/2+l)*f9(r)<
f((l+r)/2+l,r);
iff9(l)*f9(r)>
9.2结果
X=4.71239、8.24668、17.6715、14.1372
第10题
10.1程序
n=3;
%节点个数
Xj=0:
1/n:
1;
y=sin(pi*Xj);
n+1
f(i,1)=y(i);
forj=2:
fori=1:
n-j+2;
dx=((j-1)/n);
f(i,j)=(f(i+1,j-1)-f(i,j-1))/dx;
a(i)=f(1,i);
x=0:
0.001:
1/0.001+1;
y1(i)=sin(pi*x(i));
y2(i)=f10(a,Xj,x(i));
figure(10);
10.2结果
10.3分析
有图像可知插值函数的值已经很接近原函数的值了。
第11题
11.1程序
请输入n:
)%n代表节点
Xj=-5:
5;
Yj=1./(1.+Xj.^2);
x=-5:
0.01:
10/0.01+1;
y1(i)=1/(1+x(i)^2);
y2(i)=f(Yj,Xj,x(i));
figure(11);
plot(x,y2,,'
functiony=f(Yj,Xj,x)
y=0;
length(Yj)
l=1;
length(Xj)
ifi==j
continue;
end
l=l*(x-Xj(j))/(Xj(i)-Xj(j));
y=y+Yj(i)*l;
11.2结果
从左往右n依次取1、2、3、4
11.3分析
随着n的不断增加,插值越来越接近真实值。
第12题
12.1程序
)%n=501002005001000
2*pi/n:
2*pi;
X=x(i);
Y(i)=exp(3*X)*cos(pi*X);
Y
(1)=1;
X=(x(i)+x(i+1))/2;
Y2(i)=exp(3*X)*cos(pi*X);
T=(x(n+1)-x
(1))/(2*n)*(Y
(1)+2*sum(Y(2:
n))+Y(n+1));
S=(x(n+1)-x
(1))/(6*n)*(Y
(1)+4*sum(Y2)+2*sum(Y(2:
13.2结果
50
100
200
500
1000
T
3.512534119493697e+07
3.520489119998542e+07
3.522553497836321e+07
3.523136936592422e+07
3.523220478232005e+07
S
3.523140786833490e+07
3.523241623782247e+07
3.523247916807637e+07
3.523248325445201e+07
3.523248335509114e+07
第13题
13.1程序
G2=f(1/sqrt(3))+f(-1/sqrt(3));
G3=0.555555556*f(-0.7745966692)+0.555555556*f(+0.7745966692)+0.88888888889*f(0);
G5=0.2369268851*f(-0.9061798459)+0.2369268851*f(+0.9061798459)+...0.4786286705*f(-0.5384693101)+0.4786286705*f(+0.5384693101)+...0.5688888889*f(0);
f.m文件1
functiony=f(x)
y=x^2/sqrt(1-x^2);
f.m文件2
x=pi/2*(x+1)/2;
y=sin(x)/x;
13.2结果
(1)y=x^2/sqrt(1-x^2)
(2)y=sin(x)/x
G2
0.816496*********
1.74487172393660
G3
1.05409255403515
1.74531021550858
G5
1.74530859901330
第14题
14.1程序
U0=2;
%初值
step=0.2;
%step表示步长
U1
(1)=U0;
Len=1;
floor(Len/step)
U1(i+1)=U1(i)+step*f14(U1(i),(i-1)*step);
%改进的Euler法
U2
(1)=U0;
U2(i+1)=U2(i)+step/2*(f14(U2(i),(i-1)*step)+f14(U1(i+1),(i+1-1)*step));
%Runge_KUtta
U3
(1)=U0;
t=(i-1)