Matlab第2次作业数值.docx
- 文档编号:7606193
- 上传时间:2023-01-25
- 格式:DOCX
- 页数:23
- 大小:130.49KB
Matlab第2次作业数值.docx
《Matlab第2次作业数值.docx》由会员分享,可在线阅读,更多相关《Matlab第2次作业数值.docx(23页珍藏版)》请在冰豆网上搜索。
Matlab第2次作业数值
1•设a=(1,2,3),b=(2,4,3),分别计算a./b,a.\b,a/b,a\b,分析结果的意义。
解:
>>a=[1,2,3];
>>b=[2,4,3];
>>a./b【意义】矩阵内的元素对应的作-运算
b
ans=
0.50000.50001.0000
>>a.\b【意义】矩阵内的元素对应的作-运算
b
ans=
221
>>a/b【意义】矩阵整体xb=a,求x;
ans=
0.6552
>>a\b
【意义】矩阵整体ax=b,求x;
ans=
000
000
0.66671.33331.0000
2•用矩阵除法解下列线性方程组,并判断解的意义
411%9
1326X22;
153x31
解:
>>a=[41-1;32-6;1-53];
>>b=[9;-2;1];
>>a\b
ans=
2.3830
1.4894
2.0213
【意义】最小二乘法求解ax=b的近似解
Xi
X2
解:
>>a=[41;32;1-5];
>>b=[1;1;1];
>>a\bans=
0.3311
-0.1219
【意义】最小二乘法求解ax=b的近似解
2
11
1
%
X2
1
3
1
21
1
x3
2
1
12
1
3
x4
解:
>>
a=[21-11;1
21
-1;1
121
>>
b=[1;2;3];
>>a\b
ans=
1
0
1
0
【意义】最小二乘法求解欠定方程ax=b的解
3.(人口流动趋势)对城乡人口流动作年度调查,发现有一个稳定的朝向城镇流动的趋势,每年农村居民的5赠居城镇而城镇居民的1%i出,现在总人口的20%u于城镇。
加入城乡总人口保持不变,并且人口流动的这种趋势继续下去,那么
(1)一年以后住在城镇人口所占比例是什么?
两年以后呢?
十年以后呢?
(2)很多年以后呢?
(3)如果现在总人口70%u于城镇,很多年以后城镇人口所占比例是什么?
(4)计算转移矩阵的最大特征值级对应的特征向量,与问题
(2),(3)有何关系?
解:
(1)
>>a=[0.990.05;0.010.95];x0=[0.20.8]';
>>x1=a*x0,x2=aA2*x0,x10=aA10*x0
x1=
0.2380
0.7620x2=
0.2737
0.7263x10=
0.4922
0.5078
(2)
【方法一:
循环的方法】
>>x=x0;fori=1:
1000,x=a*x;
end,x
x=
0.8333
0.1667
【方法二:
累乘的方法】
>>x=aA1000*x0
x=
0.8333
0.1667
【注】若求ax=x,即(a-I)x=0的非零解,得结果如下
>>clear
a=[0.990.05:
0.010.95]:
b=*y^(2,2):
-0*0606
-0.1M1
错误原因在于,非零解是一组基础解系,不能具体确定x的值。
事实上,所求
X=Kx.(与下(4)所问对应)
(3)
>>x0=[0.70.3]';
>>x=x0;fori=1:
1000,x=a*x;end,x
x=
0.8333
0.1667
(4)>>[v,d]=eig(a)
v=
0.9806-0.7071
0.19610.7071
d=
1.00000
00.9400
>>v(:
1)./x
ans=
1.1767
1.1767
最大特征值1,[0.8333,0.1667]是对应特征值向量之一。
4.(经济预测)在某经济年度内,各经济部门的投入产出表如表3.5(单位:
亿
元)。
消耗部门
最后需求
总产值
工业
农业
第三产业
生产部门
工业
6
2
1
16
25
农业
2.25
1
0.2
1.55
5
第三产
业
3
0.2
1.8
15
20
假设某经济年度工业、农业及第三产业的最后需求均为17亿元,预测该经济年
度工业、农业及第三产业的产出(提示:
对于一个特定的经济系统而言,直接消耗矩阵和Leontief矩阵可视作不变)。
解:
>>B=[621;2.2510.2;30.21.8];x=[25520];
>>C=B/diag(x)
C=
0.24000.40000.0500
0.09000.20000.0100
0.12000.04000.0900
>>A=eye(3,3)-C
0.7600-0.4000-0.0500
-0.09000.8000-0.0100
-0.1200-0.04000.9100>>D=[171717]';x=A\Dx=
37.5696
25.7862
24.7690
所以工业37.5696,农业25.7862,第三产业24.7690.
5.
求下列矩阵的行列式、逆、特征值和特征向量。
5
765
41
1
7
1087
1
32
6
2
6
8109
15
3
5
7910
5
6
1
5
6
4
n阶方阵
1
5■
**°
n分别为5,50和500.
11
+<
6
1
5
解:
(1)
>>a=[41-1;32-6;1-53];det(a),inv(a),[v,d]=eig(a)
ans=
-94ans=
0.2553-0.02130.0426
0.1596-0.1383-0.2234
0.1809-0.2234-0.0532v=
0.0185-0.9009-0.3066
-0.7693-0.1240-0.7248
-0.6386-0.41580.6170
d=
-3.052700
03.67600
008.3766
(2)
>>a=[5765;71087;68109;57910];det(a),inv(a),[v,d]=eig(a)ans=
1ans=
68.0000-41.0000-17.000010.0000
-41.000025.000010.0000-6.0000
-17.000010.00005.0000-3.0000
d=
0.0102
00
0
0
0.84310
0
0
03.8581
0
0
0030.2887
(3)
【方法一:
循环的方法】
>>n=50;
>>a=zeros(n,n);
>>fori=1:
n,a(i,i)=5;end
>>fori=1:
(n-1),a(i,i+1)=6;end
>>fori=1:
(n-1),a(i+1,i)=1;end
>>a
【方法二:
课上所授稀疏矩阵的方法】通用稀疏矩阵:
>>n=5;>>a1=sparse(1:
n,1:
n,5*ones(1,n),n,n);
>>a2=sparse(2:
n,1:
n-1,ones(1,n-1),n,n);
>>a3=sparse(1:
n-1,2:
n,6*ones(1,n-1),n,n);
>>a=a1+a2+a3;
对角带稀疏矩阵:
>>n=5;
>>a=spdiags([ones(n,1),5*ones(n,1),6*ones(n,1)],[-1,0,1],n,n)>>a
a=
(1,1)5
(2,1)1
(1,2)6
(2,2)5
(3,2)1
(2,3)6
(3,3)5
(4,3)1
(3,4)6
(4,4)5
(5,4)1
(4,5)6
(5,5)5
后同:
>>det(a)
ans=
665
>>inv(a)
ans=
0.3173-0.58651.0286-1.62411.9489
-0.09770.4887-0.85711.3534-1.6241
0.0286-0.14290.5429-0.85711.0286
-0.0075
0.0376
-0.1429
0.4887
-0.5865
0.0015
-0.0075
0.0286
-0.0977
0.3173
>>[v,d]=eig(a)
v=
-0.7843
-0.7843
-0.9237
0.9860
-0.9237
0.5546
-0.5546
-0.3771
-0.0000
0.3771
-0.2614
-0.2614
0.0000
-0.1643
-0.0000
0.0924
-0.0924
0.0628
-0.0000
-0.0628
-0.0218
-0.0218
0.0257
0.0274
0.0257
0.7574
0
0
0
0
0
9.2426
0
0
0
0
0
7.4495
0
0
0
0
05.0000
0
0
0
0
02.5505
将n=50和500带入,类似可求得
6.
(1)判断第5题各小题是否可以相似对角化,如果可以,求出对角矩阵和对应的相似变换矩阵。
解:
若n阶矩阵含有n个不等的特征值,那么他一定可以相似对角化,且此n个特征向量组成的矩阵即相似变换矩阵。
故第
(1)题、第
(2)题、第(3)题的n=5时,可以相似对角化。
对角矩阵就是特征值矩阵,相似变换矩阵就是特征向量矩阵。
(2)判断第5题各小题是否为正定矩阵。
解:
【XX】matlab里如何判断矩阵为正定:
[Dp]=chol(A),如果A正定,返回的p=0,如果不正定,则返回一个正的p,p—1为A中正定子矩阵的阶次,即D为p—1为D的阶次。
>>a=[41-1;32-6;1-53];
>>[Dp]=chol(a)
D=
2.00000.5000
01.3229
p=
3
P工0,故矩阵不正定。
>>a=[5765;71087;68109;57910];[Dp]=chol(a)
D=
2.23613.13052.68332.2361
00.4472-0.89440
001.41422.1213
0000.7071
0
P=0,故矩阵正定>>n=5;
>>a=spdiags([ones(n,1),5*ones(n,1),6*ones(n,1)],[-1,0,1],n,n);
>>[Dp]=chol(a)
D=
(1.1)2.2361
(1.2)2.6833
p=
2
Pm0,故矩阵不正定。
7•求下列向量组的秩和它的一个最大线性无关组,并将其余向量用该最大无关组线性表示。
1=7,-6,-7,0.
1=4,-3,1,32=2,-13,5,3=1,-1,-1,1,4=3,-2,3,4,
解:
>>a=[4-313;2-135;1-1-1-1;3-234;7-6-70];
>>rank(a)ans=
3
>>rank(a(1:
3,:
))ans=
>>rank(a([124],:
))ans=
>>b=a([124],:
)';c=a([35],:
)';>>b\cans=
0.50005.0000
-0.50001.0000
-0.0000-5.0000
[a1a2a4]构成了一个极大线性无关组。
由计算结果可以知道,
a3=0.5a-0.5a2,a5=5a1+a2-5a4.
8•假定某天的气温变化记录如下表,试用最小二乘方法找出这一天的气温变化规律.
时刻助)
0
1
J
E
1
5
6
7
s
9
10
11
12
温度
w
15°
14°
uc
14°
U°
16°
23°
25°
2S°
时刻
13
14
15
16
17
18
19
20
21
)2
23
24
31fl
31°
29°
27°
2中
24°
■
2(f
17°
16°
解:
>>a=0:
24;
>>b=[15141414141516182022232528313231292725242220181716];
>>a仁polyfit(a,b,1);
>>a2=polyfit(a,b,2);
>>a3=polyfit(a,b,3);
>>a4=polyfit(a,b,4);
>>a5=polyfit(a,b,5);
>>y1=polyval(a1,a);
>>plot(a,b,'r*',a,y1,'b-');
>>y2=polyval(a2,a);
>>y3=polyval(a3,a);
>>y4=polyval(a4,a);
>>y5=polyval(a5,a);
>>plot(a,b,'r*',a,y2,'b-');
>>plot(a,b,'r*',a,y3,'b-');
>>plot(a,b,'r*',a,y4,'b-');
>>plot(a,b,'r*',a,y5,'b-');
>>a5
a5=
0.0001-0.00370.04500.0157-0.830415.0539
拟合曲线如下:
9.一矿脉有13个相邻样本点,人为地设定一原点,现测得各样本点对原点的距离x,与样本点处某种金属含量y的一组数据如下,画出散点图观察二者的关
系,
试建立合适的模型,
如二次曲线
、双曲线、对数曲线等。
X
V
rJ.
4
57
8
10
y
I06J2
109J0
109耶
109.50
11D.00
109,93
110.19
X
JI
14
15
15
18
19
V
i|T
110.59
11060
110.90
110.76
111.00
111.20
二次曲线(y=ax2+bx+c):
>>x=[23457810111415151819];
>>y=[106.42109.2109.58109.5110109.93110.49110.59110.6110.9110.76111111.2];
>>a=polyfit(x,y,2);
>>y1=polyval(a,x);
>>plot(x,y1,'b-',x,y,'r*')
>>c1=y1-y;
>>dot(c1,c1)ans=
4.1821
112
oo
双曲线(y=-+?
?
:
>>a=[1./x',ones(13,1)]
>>ab=a\y'
ab=
-9.0300
111.4405>>y2=ab
(1)./x+ab
(2);
>>plot(x,y2,'b-',x,y,'r*')c2=y2-y;dot(c2,c2)
ans=
对数曲线(y=alogx+b):
>>a=[log(x'),ones(13,1)]>>ab=a\y'
ab=
1.5663
106.7113
>>y3=ab
(1).*log(x)+ab
(2);
>>plot(x,y3,'b-',x,y,'r*')c3=y3-y;
dot(c3,c3)
ans=
双曲线拟合效果最佳
10.求下列积分的数值解
2
(1)e2xcos3(x)dx
0
(2)(1xy2)dydx,D为x2y22x
D
解:
(1)
fun.m
functiony=fun(t)y=exp(2*t).*((cos(t)).A3);
>>d=pi/1000;
>>t=0:
d:
2*pi;
>>nt=length(t);
>>y=fun(t);
>>sc=cumsum(y)*d;
>>scf=sc(nt)
scf=
9.7505e+04
>>z=trapz(y)*d
z=
9.7054e+04
>>quad('fun',0,2*pi)
ans=
9.7054e+04
矩形法误差较大,梯形和simpson法结果比较精确
(2)
>>fun=@(x,y)1+x+y42;
>>ymax=@(x)(2.*x-x.A2).A0.5;
>>ymin=@(x)-((2.*x-x.A2).A0.5);
>>quad2d(fun,0,2,ymin,ymax)
11.用积分法计算下列椭园的周长
解:
所求周长,利用椭圆的参数方程,等价于求积分
2
■■-45cos2d
0
>>fun=@(x)(4+5*(cos(x)).A2).A0.5;
>>quad(fun,0,2*pi)ans=
15.8654
12.试求下列积分,出现什么问题?
分析原因,设法求出正确的解
1
Ix0.2cos(x)dx
1
解:
>>fun=@(x)(x.A0.2).*cos(x);
>>quad(fun,-1,1)
ans=
1.2390+0.4026i
【问题】所得答案为虚数,明显错误。
【错因】在matlab里面,对于区间(-1,0),xA0.2结果是虚数(数学上是实数才对!
!
具体原因不知啊,求老师课上解答)
【正解】因为x0.2cos(x)为奇函数,积分区间(-1,1)关于原点对称,所以
1=0O
k
13考虑积分1sin(x)dx
0
2k试分别用trapz(取步长h=0.1或n),quad
和quadl求解1(8)和I(32)
(1)I(8):
。
发现什么问题?
先求I(8)
>>fun=@(x)abs(sin(x));
>>d=0.1;
>>t=0:
d:
8*pi;
>>y=fun(t);
>>z=trapz(y)*d
z=
15.9981
z=trapz(y)*d
z=
1.2311e-14
>>quad(fun,0,8*pi)
ans=
16.0000
d=pi;t=0:
d:
8*pi;
>>quadl(fun,0,8*pi)
ans=
y=fun(t);16.0000
quad、quadl两种方法所求结果均正确。
在间距为0.1时,trapz所得结果也比较精确。
但间距为pi时,由于所取步长过大,结果错误。
这是因为运用NewtonCotes公式时,误差太大了!
而应取远小于积分长度的间距。
*14.用正交变换化下列二次型为标准形
222f(x1,X2,X3)=X1-4X1X2+4X1x3-2x2+8x2X3-2x3
>>a=[1-22;-2-24;24-2];
>>[v,d]=eig(a)
v=
0.33330.9339-0.1293
0.6667-0.3304-0.6681
-0.66670.1365-0.7327
d=
-7.000000
02.00000
002.0000
>>v'*v
ans=
1.00000.00000.0000
0.00001.00000
0.000001.0000
标准型即为f=-7y1A2+2y2A2+2y3A2
*15.下图是连接三个电压已知终端的电路网,求
由基尔霍夫方程:
1311
?
?
?
?
?
2=10
12…4…3''
11135
?
s+?
?
?
?
=-
3…5…15…3
1471
?
?
?
?
+?
?
=04…60…5''
>>a=[13/12-1/4-1/3;1/31/5-13/15;1/4-47/601/5];
>>b=[10;-5/3;0];
>>a\b
ans=
13.3453
6.4401
8.5420
即a=13.3453V,b=6.4401V,c=8.5420V
123
456
*16.就矩阵A=验证下列性质
780
(i)设1,2,
n为n阶方阵A的特征值,则
i1
aH(A的迹),
i1
n
i
(1)nA
i1
(ii)设f(x)为A的特征多项式,则f(A)=0解:
(i)
>>A=[123;456;780];
>>trace(A)
ans=
6
>>[v,d]=eig(A)
v=
-0.2998-0.7471-0.2763
-0.70750.6582-0.3884
-0.6400-0.09310.8791
d=
12.122900
0-0.38840
00-5.7345
>>trace(d)
ans=
6.0000>>det(A)
ans=
27.0000>>det(d)ans=
27.0000(ii)
>>root1=[d(1,1)d(2,2)d(3,3)];
>>p=poly(root1)
P=
1.000-6.0000-72.0000-27.0000>>polyvalm(p,A)ans=
1.0e-11*
-0.0373-0.0462-0.0380
-0.0853-0.1133-0.0860
-0.0853-0.1137-0.0625
可以发现,结果约等于0.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 作业 数值