MATLAB教程 R2014a 答案 全 张志涌-matlab2014答案
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
目录
第一章 (1)
第二章 (5)
第三章 (12)
第四章 (32)
第五章 (47)
第六章 (54)
补充题欧拉法,龙格库塔法解方程,黑板上的题 (57)
第一章
1.创建表达式
%可以用syms先符号运算再带入值
x=1;
y=2;
z=(sqrt(4*x^2+1)+0.5457*exp(-0.75*x^2-3.75*y^2-1.5*x))/(2*sin(3*y)-1)
z =
-1.4345
2.计算复数
x=(-1+sqrt(-5))/4;
y=x+8+10j
y =
7.7500 +10.5590i
3.help命令学三维曲线
x=-5:0.1:5;
y=x;
[X,Y]=meshgrid(x,y);
Z=(sin(sqrt(X.^2+Y.^2)))./(sqrt(X.^2+Y.^2));
subplot(221);
surf(X,Y,Z);
colormap(cool);
subplot(222);
plot3(X,Y,Z,'linewidth',4); %绘制三维曲线,也可以随意给定一个三维曲线的函数。
如果画这个曲面,那么将绘出一族三维曲线
grid on;
subplot(223);
meshz(X,Y,Z); %地毯绘图
subplot(224);
meshc(X,Y,Z); %等高线绘图
4.peaks等高线(更改原函数)
subplot(221);
contour(peaks1,20);
subplot(222);
contour3(peaks1,10); %可以定义等高线条数
subplot(223);
contourf(peaks1,10);
subplot(224);
peaks1;
z = 3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) ...
- 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...
- 1/3*exp(-(x+1).^2 - y.^2)
5. LOGO绘制membrane
logo
第一章书后习题
1.合法性
不合法合法不合法不合法合法
2.运行命令及探讨
a=sqrt(2)
a =
1.4142。
是一个近似。
可通过改变format进行位数显示调整。
例如:
format long;
a=sqrt(2)
format short;
a =
1.414213562373095
或可使用digits任意指定输出位数。
例如:
digits(50);
a=sqrt(2);
vpa(a)
ans =
1.4142135623730950488016887242096980785696718753769
常见情况下毋需太高精度。
3.运行结果讨论
format long;
w1=a^(2/3)
w2=a^2^(1/3)
w3=(a^(1/3))^2
w1 =
1.259921049894873
w2 =
1.259921049894873
w3 =
1.259921049894873
测试结果为相同,说明MATLAB程序执行时经过的过程相同。
4.clear clf clc
clear 为从内存中清除变量和函数
clf 为清除figure中的已绘图形以及子图形
clc 为清除命令行窗口
5.产生二维数组
显然第一第二个方法可以实现。
例如:
s=[1 2 3;4 5 6;7 8 9]
s =
1 2 3
4 5 6
7 8 9
即是一个简便的键入矩阵的方法。
第二章
1 数据类型
class(3/7+0.1)
class(sym(3/7+0.1))
class(vpa(sym(3/7+0.1),4))
class(vpa(sym(3/7+0.1)))
ans =
double
ans =
sym
ans =
sym
ans =
sym
2 哪些精准?
a1=sin(sym(pi/4)+exp(sym(0.7)+sym(pi/3)));
a2=sin(sym(pi/4)+exp(sym(0.7))*exp(sym(pi/3)));
a3=sin(sym('pi/4')+exp(sym('0.7'))*exp(sym('pi/3'))); a4=sin(sym('pi/4')+exp(sym('0.7+pi/3')));
a5=sin(sym(pi/4)+exp(sym(0.7+pi/3)));
a6=sin(sym(pi/4)+sym(exp(0.7+pi/3)));
a7=sin(sym(pi/4+exp(0.7+pi/3)));
a8=sym(sin(pi/4+exp(0.7+pi/3)));
digits(64);
vpa(a2-a1)
vpa(a3-a1)
vpa(a4-a1) %为精确值
vpa(a5-a1)
vpa(a6-a1)
vpa(a7-a1)
vpa(a8-a1)
ans =
8.772689107613377606024459313047548287536202098197290121158158175e-72 ans =
8.772689107613377606024459313047548287536202098197290121158158175e-72 ans =
0.0
ans =
-0.00000000000000088748227169595846195226372540142491282548756502081529 37300697045
ans =
-0.00000000000000148912212817656334175571371627278077803022761502222373 5634526288
ans =
-0.00000000000000151885559392782263589708294774441179495071438346616836 4259064934
ans =
-0.00000000000000151859755909122793880734918235619076228065004813152159 311456667
可以看到,除了a4为精确,其余均存在很小的误差。
其中a2与a3的误差较小,小于eps 精度,故可认为为精确的。
3 独立自由变量
a1=sym('sin(w*t)') ;
a2=sym('a*exp(-X)' );
a3=sym('z*exp(j*th)');
symvar(a1,1)
symvar(a2,1)
symvar(a3,1)
ans =
w
ans =
a
ans =
z
6 符号解
syms x k;
f1=x.^k;
s1=symsum(f1,k,0,inf);
s2=subs(f1,x,(-1/3));
s3=subs(f1,x,(1/pi));
s4=subs(f1,x,3);
symsum(s2,k,0,inf)
double(symsum(s3,k,0,inf))
symsum(s4,k,0,inf)
ans =
3/4
ans =
1.4669
ans =
Inf
7 限定性假设
reset(symengine);
syms k;
syms x positive;
f1=(2/(2*k+1))*((x-1)/(x+1))^(2*k+1);
f1_s=symsum(f1,k,0,inf);
simplify(f1_s,'steps',27,'IgnoreAnalyticConstraints',true)
ans =
log(x)
8 符号计算
syms t;
yt=abs(sin(t));
dydt=diff(yt,t)
dydt0=limit(dydt,t,0,'left')
dydtpi=subs(dydt,t,(pi/2))
dydt =
sign(sin(t))*cos(t)
dydt0 =
-1
dydtpi =
9 积分值
syms x;
fx=exp(-abs(x))*abs(sin(x))
fxint=int(fx,-5*pi,1.7*pi);
vpa(fxint,64)
fx =
abs(sin(x))*exp(-x)
ans =
3617514.635647088707100018393465500554242735057835123431773680704 10二重积分
syms x y;
fxy=x^2+y^2;
int(int(fxy,y,1,x^2),x,1,2)
ans =
1006/105
11 绘出曲线
syms t x;
fx=int((sin(t)./t),t,0,x);
ezplot(fx)
fx4=subs(fx,x,4.5)
fx4 =
12 积分表达式
syms x;
syms n positive;
yn=int((sin(x)).^n,x,0,pi/2)
yn3=subs(yn,n,1/3);
vpa(yn3,32)
yn =
beta(1/2, n/2 + 1/2)/2
ans =
1.2935547796148952674767575125656
13 序列卷积
syms a b n;
syms k positive;
xk=a.^k;
hk=b.^k;
kn=subs(xk,k,k-n)*subs(hk,k,n);
yk=symsum(kn,n,0,k)
yk =
piecewise([a == b and b ~= 0, b^k*(k + 1)], [a ~= b or b == 0, (a*a^k - b*b^k)/(a - b)])
所以答案为a*a^k - b*b^k)/(a - b)
20求解solve
reset(symengine)
syms x y;
s=solve('x^2+y^2-1','x*y-2','x','y')
s.x
s.y
s =
x: [4x1 sym]
y: [4x1 sym]
ans =
((15^(1/2)*i)/2 + 1/2)^(1/2)/2 - ((15^(1/2)*i)/2 + 1/2)^(3/2)/2 - ((15^(1/2)*i)/2 + 1/2)^(1/2)/2 + ((15^(1/2)*i)/2 + 1/2)^(3/2)/2 (1/2 - (15^(1/2)*i)/2)^(1/2)/2 - (1/2 - (15^(1/2)*i)/2)^(3/2)/2 - (1/2 - (15^(1/2)*i)/2)^(1/2)/2 + (1/2 - (15^(1/2)*i)/2)^(3/2)/2 ans =
((15^(1/2)*i)/2 + 1/2)^(1/2)
-((15^(1/2)*i)/2 + 1/2)^(1/2)
(1/2 - (15^(1/2)*i)/2)^(1/2)
-(1/2 - (15^(1/2)*i)/2)^(1/2)
23 求通解
clear all;
yso=simplify(dsolve('Dy*y*0.1+0.3*x=0','x'))
yso =
(- 3*x^2 + 2*C3)^(1/2)
-(- 3*x^2 + 2*C3)^(1/2)
%此题存疑
hold on;clear all;
reset(symengine);
syms x;
y1=(- 3*x^2 + 2*1)^(1/2);
y2=-(- 3*x^2 + 2*1)^(1/2);
h1=ezplot(y1,x,[-2 2 -2 2],1);
h2=ezplot(y2,x,[-2 2 -2 2],1);
grid on;title('');warning off;axis([-2 2 -2 2]);
set(h1,'color','r','linewidth',2);
set(h2,'color','r','linewidth',2);
xlabel('Y');ylabel('x');
%对于第二章存在问题的习题的探讨2.23
clear all;
syms x;
yso=simplify(dsolve('Dy*y*0.1+0.3*x=0','x'))
%此题存疑
hold on;clear all;
reset(symengine);
syms x;
y1=(- 3*x^2 + 2*1)^(1/2);
y2=-(- 3*x^2 + 2*1)^(1/2);
h1=ezplot(y1,x,[-2 2 -2 2],1);
h2=ezplot(y2,x,[-2 2 -2 2],1);
grid on;title('');warning off;axis([-2 2 -2 2]);
set(h1,'color','r','linewidth',2);
set(h2,'color','r','linewidth',2);
xlabel('Y');ylabel('x');
yso =
(- 3*x^2 + 2*C3)^(1/2)
-(- 3*x^2 + 2*C3)^(1/2)
%以上方法可以绘出正常的横坐标为y纵坐标为x的图像,但发现在y=0处x延伸至正负无穷。
h1=ezplot(y1,[-2 2 -2 2],1);
h2=ezplot(y2,[-2 2 -2 2],1);
%以上方法绘出的图像存在一个空隙,且默认为y-x图像。
reset(symengine);
syms x y S;
S = dsolve('Dy*y/5+x/4=0','x');
ezplot(subs(y^2-(S(1))^2, 'C3', 1),[-2,2 -2,2],2); grid on;
%用椭圆方程绘图不产生间隙
24 一阶微分方程
syms a b;
ys=dsolve('Dy-a*x^2-b*x=0','y(0)=2',x)
ys =
(x^2*(3*b + 2*a*x))/6 + 2
25 边值问题
fs=dsolve('Df-3*f=4*g,Dg+4*f=3*g','f(0)=0,g(0)=1') fs =
g: [1x1 sym]
f: [1x1 sym]
fs.g
fs.f
ans =
cos(4*t)*exp(3*t)
ans =
sin(4*t)*exp(3*t)
第三章
3.行下标列下标
rng('default');
A=rand(3,5);
L=A>0.5
L =
1 1 0 1 1
1 1 1 0 0
0 0 1 1 1
[a,b]=find(L==1)
IND=sub2ind(size(A),a,b)
IND =
1
2
4
5
8
9
10
12
13
15
4.循环运算、数组运算
t=0:0.1:10;
N=length(t);
y1=zeros(size(t));
for k=1:N
y1(k)=1-exp(-0.5*t(k))*cos(2*t(k)); end
plot(t,y1);
xlabel('t');
ylabel('y1');
y2=1-exp(-0.5*t).*cos(2*t); plot(t,y2);
xlabel('t');
ylabel('y2');
5.回答问题
clear all;
A=magic(3);B=rand(3);
A*B
B*A
ans =
5.4072 11.5771 3.0037
6.3884 10.3215 4.9680
2.7058 7.5337 4.8496
ans =
2.5916
3.8303 5.2097
3.4833 5.6313 3.6800
10.9646 9.0086 12.3554
相同,对于矩阵而言对位相乘无差异
不相同,点乘与矩阵乘法进行的不是同一种运算。
不相同,左乘右乘运算不同。
相同,A左点除B等同于B右点除A,均是对位计算。
不相同,左除右除运算亦不相同。
A*A\B-B
A*(A\B)-B
A*(A*inv(B))-B
ans =
-0.0562 -0.6902 -0.0436
-0.1051 -0.3282 -0.4311
-0.8011 -0.9350 -0.3763
ans =
1.0e-15 *
0 -0.1110 -0.0278
0 0 0
0 0.1110 0
ans =
-80.2971 65.0383 107.2212
-8.0299 91.2626 70.5679
-66.2535 153.4898 66.4342
不相同。
第二个更接近0。
具体原理需要参考线性代数书……有点忘了。
A\eye(3)
eye(3)/A
ans =
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
ans =
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
相同。
因为对于对角阵,,二者均可化为同一形式。
6.结果不同
A=[1 2; 3 4];
B1=A.^(0.5)
B2=0.5.^A
B3=A^(0.5)
B4=0.5^A
B1 =
1.0000 1.4142
1.7321
2.0000
B2 =
0.5000 0.2500
0.1250 0.0625
B3 =
0.5537 + 0.4644i 0.8070 - 0.2124i
1.2104 - 0.3186i 1.7641 + 0.1458i
B4 =
0.9910 -0.4422
-0.6634 0.3276
A1=B1.*B1
A3=B3*B3
norm(A1-A3,'fro')
A1 =
1.0000
2.0000
3.0000
4.0000
A3 =
1.0000 + 0.0000i
2.0000 + 0.0000i
3.0000 - 0.0000i
4.0000 + 0.0000i ans =
1.2831e-15
可见误差在eps量级,可以认为相等。
7.绘出图形
x=-3*pi:pi/15:3*pi;
y=x;
[X,Y]=meshgrid(x,y);
warning off;
Z=sin(X).*sin(Y)./X./Y;
共有10个非数数据。
surf(X,Y,Z)
shading interp
x=-3*pi:pi/15:3*pi;
Lx=(x==0);
xx=x+Lx*realmin;
y=xx;
[X,Y]=meshgrid(xx,y);
warning off;
Z=sin(X).*sin(Y)./X./Y;
surf(X,Y,Z)
shading interp
即消除零点处的断点即可
8.两种思路
%第二种思路
function z=zpoly_z(x,y)
if x+y<=-1
z=0.546*exp(-0.75*y.^2-3.75*x.^2+1.5*x); elseif x+y>-1 & x+y<=1
z=0.758*exp(-y.^2-6*x.^2);
else
z=0.546*exp(-0.75*y.^2-3.75*x.^2-1.5*x); end
x=-1.5:0.1:1.5;
y=-3:0.1:3;
[X,Y]=meshgrid(x,y);
Z=zpoly_z(X,Y);
surf(X,Y,Z);
%第一种思路
x=-1.5:0.1:1.5;
y=-3:0.2:3;
LX=length(x);
LY=length(y);
for ii=1:LX
for jj=1:LY
if x(ii)+y(jj)<=-1
z=0.546*exp(-0.75*y.^2-3.75*x.^2+1.5*x); elseif x(ii)+y(jj)>-1 & x(ii)+y(jj)<=1
z=0.758*exp(-y.^2-6*x.^2);
else
z=0.546*exp(-0.75*y.^2-3.75*x.^2-1.5*x); end
end
end
[X,Y]=meshgrid(x,y);
Z=zpoly_z(X,Y);
surf(X,Y,Z);
%其实for循环完全无意义……
9.矩阵计算
%第一问老师取消
rng default
A=randn(50,70)+1i*randn(50,70); B=randn(70,60)+1i*randn(70,60); C=randn(50,60)+1i*randn(50,60); D=randn(60,1)+1i*randn(60,1);
G=(A*B-C)*D
Gr=real(G),70,70
Gi=imag(G)
Gn=norm(G,2)
G =
1.0e+02 *
-0.1776 + 1.9914i
0.6088 + 0.3316i
-0.1340 - 0.8615i
0.0752 - 0.0759i
-0.1171 - 1.8169i
0.2005 - 1.4540i
-1.4501 + 0.1897i
0.6445 + 0.1657i
-1.0651 + 0.1191i
0.3301 - 0.0450i
-1.4338 + 0.8707i
-0.9491 + 1.4840i
1.1314 + 1.2751i
-0.5158 - 0.0725i
-0.2746 + 0.2518i
-1.0279 - 0.8409i
-1.1161 - 2.3362i 0.1346 + 1.3500i 0.4220 - 1.2839i 0.2650 - 0.2849i -1.0212 + 0.5374i 0.0563 + 0.4151i -1.9074 - 0.2448i
0.1645 + 1.2071i
1.1870 + 0.0085i 1.2304 + 0.6672i 0.3303 - 1.6027i -0.5728 - 0.5519i 0.3738 + 0.2863i -0.6682 - 0.7565i 1.6063 + 1.2886i 0.6994 - 1.3377i 0.6523 + 0.0318i -0.2143 -
2.8209i 1.7026 - 0.1371i 0.9285 + 1.5852i -0.7550 - 0.2427i -1.3879 - 1.8978i -0.5266 - 0.8334i -0.0849 + 0.1680i 1.1590 + 0.2109i -1.8938 + 0.6709i 0.3406 - 1.8211i -1.0916 - 1.8076i
0.2062 - 1.4363i
1.3679 + 0.2061i -0.4541 + 0.8056i 1.3574 + 0.8773i -0.1071 + 0.0948i 0.1042 +
2.2812i Gr =
-17.7553
60.8848
-13.4003
7.5175
-11.7073
20.0458
-145.0055
64.4517
-106.5069
33.0077
-143.3779
-94.9055
113.1368
-51.5804
-27.4560
-102.7914
-111.6150
13.4596
42.2009
26.5006
-102.1225
5.6295
-190.7388 16.4525 118.6963 123.0361 33.0336 -57.2817 37.3849 -66.8175 160.6261 69.9436 65.2278 -21.4319 170.2597 92.8549 -75.5045 -138.7923 -52.6574 -8.4902 115.9030 -189.3844 34.0593 -109.1584 20.6169 136.7896 -45.4089 135.7386 -10.7050 10.4240 Gi =
199.1404 33.1590 -86.1452 -7.5887 -181.6856 -145.4039 18.9686 16.5731 11.9053 -4.5021 87.0651 148.4022 127.5072 -7.2483 25.1791 -84.0887 -233.6194 135.0018 -128.3931 -28.4923 53.7385 41.5139 -24.4788 120.7113 0.8532 66.7238 -160.2738 -55.1871
28.6287
-75.6522
128.8596
-133.7671
3.1772
-282.0866
-13.7111
158.5203
-24.2673
-189.7767
-83.3384
16.7992
21.0869
67.0898
-182.1134
-180.7631
-143.6344
20.6149
80.5622
87.7339
9.4764
228.1237
Gn =
1.0253e+03
y2=1-exp(-0.5*t).*cos(2*t); plot(t,y2);
xlabel('t');
ylabel('y2');
5.回答问题
clear all;
A=magic(3);B=rand(3);
A*B
B*A
ans =
5.4072 11.5771 3.0037
6.3884 10.3215 4.9680
2.7058 7.5337 4.8496
ans =
2.5916
3.8303 5.2097
3.4833 5.6313 3.6800
10.9646 9.0086 12.3554
相同,对于矩阵而言对位相乘无差异
不相同,点乘与矩阵乘法进行的不是同一种运算。
不相同,左乘右乘运算不同。
相同,A左点除B等同于B右点除A,均是对位计算。
不相同,左除右除运算亦不相同。
A*A\B-B
A*(A\B)-B
A*(A*inv(B))-B
ans =
-0.0562 -0.6902 -0.0436
-0.1051 -0.3282 -0.4311
-0.8011 -0.9350 -0.3763
ans =
1.0e-15 *
0 -0.1110 -0.0278
0 0 0
0 0.1110 0
ans =
-80.2971 65.0383 107.2212
-8.0299 91.2626 70.5679
-66.2535 153.4898 66.4342
不相同。
第二个更接近0。
具体原理需要参考线性代数书……有点忘了。
A\eye(3)
eye(3)/A
ans =
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
ans =
0.1472 -0.1444 0.0639
-0.0611 0.0222 0.1056
-0.0194 0.1889 -0.1028
相同。
因为对于对角阵,,二者均可化为同一形式。
6.结果不同
A=[1 2; 3 4];
B1=A.^(0.5)
B2=0.5.^A
B3=A^(0.5)
B4=0.5^A
B1 =
1.0000 1.4142
1.7321
2.0000
B2 =
0.5000 0.2500
0.1250 0.0625
B3 =
0.5537 + 0.4644i 0.8070 - 0.2124i
1.2104 - 0.3186i 1.7641 + 0.1458i
B4 =
0.9910 -0.4422
-0.6634 0.3276
A1=B1.*B1
A3=B3*B3
norm(A1-A3,'fro')
A1 =
1.0000
2.0000
3.0000
4.0000
A3 =
1.0000 + 0.0000i
2.0000 + 0.0000i
3.0000 - 0.0000i
4.0000 + 0.0000i ans =
1.2831e-15
可见误差在eps量级,可以认为相等。
7.绘出图形
x=-3*pi:pi/15:3*pi;
y=x;
[X,Y]=meshgrid(x,y);
warning off;
Z=sin(X).*sin(Y)./X./Y;
共有10个非数数据。
surf(X,Y,Z)
shading interp
x=-3*pi:pi/15:3*pi;
Lx=(x==0);
xx=x+Lx*realmin;
y=xx;
[X,Y]=meshgrid(xx,y);
warning off;
Z=sin(X).*sin(Y)./X./Y;
surf(X,Y,Z)
shading interp
即消除零点处的断点即可
8.两种思路
%第二种思路
function z=zpoly_z(x,y)
if x+y<=-1
z=0.546*exp(-0.75*y.^2-3.75*x.^2+1.5*x); elseif x+y>-1 & x+y<=1
z=0.758*exp(-y.^2-6*x.^2);
else
z=0.546*exp(-0.75*y.^2-3.75*x.^2-1.5*x); end
x=-1.5:0.1:1.5;
y=-3:0.1:3;
[X,Y]=meshgrid(x,y);
Z=zpoly_z(X,Y);
surf(X,Y,Z);
%第一种思路
x=-1.5:0.1:1.5;
y=-3:0.2:3;
LX=length(x);
LY=length(y);
for ii=1:LX
for jj=1:LY
if x(ii)+y(jj)<=-1
z=0.546*exp(-0.75*y.^2-3.75*x.^2+1.5*x); elseif x(ii)+y(jj)>-1 & x(ii)+y(jj)<=1
z=0.758*exp(-y.^2-6*x.^2);
else
z=0.546*exp(-0.75*y.^2-3.75*x.^2-1.5*x); end
end
end
[X,Y]=meshgrid(x,y);
Z=zpoly_z(X,Y);
surf(X,Y,Z);
%其实for循环完全无意义……
9.矩阵计算
%第一问老师取消
rng default
A=randn(50,70)+1i*randn(50,70); B=randn(70,60)+1i*randn(70,60); C=randn(50,60)+1i*randn(50,60); D=randn(60,1)+1i*randn(60,1);
G=(A*B-C)*D
Gr=real(G),70,70
Gi=imag(G)
Gn=norm(G,2)
G =
1.0e+02 *
-0.1776 + 1.9914i
0.6088 + 0.3316i
-0.1340 - 0.8615i
0.0752 - 0.0759i
-0.1171 - 1.8169i
0.2005 - 1.4540i
-1.4501 + 0.1897i
0.6445 + 0.1657i
-1.0651 + 0.1191i
0.3301 - 0.0450i
-1.4338 + 0.8707i
-0.9491 + 1.4840i
1.1314 + 1.2751i
-0.5158 - 0.0725i
-0.2746 + 0.2518i
-1.0279 - 0.8409i
-1.1161 - 2.3362i 0.1346 + 1.3500i 0.4220 - 1.2839i 0.2650 - 0.2849i -1.0212 + 0.5374i 0.0563 + 0.4151i -1.9074 - 0.2448i
0.1645 + 1.2071i
1.1870 + 0.0085i 1.2304 + 0.6672i 0.3303 - 1.6027i -0.5728 - 0.5519i 0.3738 + 0.2863i -0.6682 - 0.7565i 1.6063 + 1.2886i 0.6994 - 1.3377i 0.6523 + 0.0318i -0.2143 -
2.8209i 1.7026 - 0.1371i 0.9285 + 1.5852i -0.7550 - 0.2427i -1.3879 - 1.8978i -0.5266 - 0.8334i -0.0849 + 0.1680i 1.1590 + 0.2109i -1.8938 + 0.6709i 0.3406 - 1.8211i -1.0916 - 1.8076i
0.2062 - 1.4363i
1.3679 + 0.2061i -0.4541 + 0.8056i 1.3574 + 0.8773i -0.1071 + 0.0948i 0.1042 +
2.2812i Gr =
-17.7553
60.8848
-13.4003
7.5175
-11.7073
20.0458
-145.0055
64.4517
-106.5069
33.0077
-143.3779
-94.9055
113.1368
-51.5804
-27.4560
-102.7914
-111.6150
13.4596
42.2009
26.5006
-102.1225
5.6295
-190.7388 16.4525 118.6963 123.0361 33.0336 -57.2817 37.3849 -66.8175 160.6261 69.9436 65.2278 -21.4319 170.2597 92.8549 -75.5045 -138.7923 -52.6574 -8.4902 115.9030 -189.3844 34.0593 -109.1584 20.6169 136.7896 -45.4089 135.7386 -10.7050 10.4240 Gi =
199.1404 33.1590 -86.1452 -7.5887 -181.6856 -145.4039 18.9686 16.5731 11.9053 -4.5021 87.0651 148.4022 127.5072 -7.2483 25.1791 -84.0887 -233.6194 135.0018 -128.3931 -28.4923 53.7385 41.5139 -24.4788 120.7113 0.8532 66.7238 -160.2738 -55.1871
28.6287
-75.6522
128.8596
-133.7671
3.1772
-282.0866
-13.7111
158.5203
-24.2673
-189.7767
-83.3384
16.7992
21.0869
67.0898
-182.1134
-180.7631
-143.6344
20.6149
80.5622
87.7339
9.4764
228.1237
Gn =
1.0253e+03
第四章
1 数值差分
d=0.001;
dxdt_diff=diff(y)/d;
dxdt_grad=gradient(y)/d;
hold on;
plot(t,y,'-r');
plot(t(1:end-1),dxdt_diff,'-g'); plot(t,dxdt_grad,'-b');
2 数值计算
d=1e-5;
t=0:d:10;
t=t+(t==0)*realmin;%去除影响fx=sin(t)./t;
y=cumtrapz(fx)*d;
plot(t,y);
yt1=find(t==4.5);
y45=y(yt1)
y45 =
1.6541
d=1e-5;
t=0:d:10;
t=t+(t==0)*realmin;%去除影响
fx=@(t)sin(t)./t;
y=integral(fx,0,10)
%plot(t,y);
%yt1=find(t==4.5);
%y45=y(yt1)
y =
1.6583
3 数值积分
d=1e-5;
t=0:d:pi;
t=t+(t==0)*realmin;%去除影响
fx=@(t)exp((sin(t)).^3);
s=integral(fx,0,pi)
s =
5.1370
符号验算
syms x;
fxx=exp((sin(x))^3);
ss=int(fxx,x,0,pi)
ss =
int(exp(sin(x)^3), x, 0, pi)
此时符号运算不如数值运算可以得出解析解
4 求数值积分
clear;
fx=@(x)exp(-abs(x)).*abs(sin(x));
format long ;
s=integral(fx,-5*pi,inf,'Abstol',1e-9)
s =
1.090331100328886
x=linspace(-5*pi,50,1e7);
dx=x(2)-x(1);
st=trapz(exp(-abs(x)).*abs(sin(x)))*dx
st =
1.090331328559288
可知,trapz必须为有限数值;同样,quadl也需为有限数值
clear;
fx=@(x)exp(-abs(x)).*abs(sin(x));
format long ;
s=quadl(fx,-5*pi,50,1e-9)
s =
1.090331336243919
5 最小值点
%首先绘出函数曲线,确定大概找最小值范围
clear;
ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t); t=-5:0.01:5;
ezplot(ft,[-5,5]);
%显然,最小值在[-2 -1]区间内。
tt=linspace(-2,-1,101);
for k=1:100
[tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1));
end
mini=min(fobj)
tmini=tmin(find(fobj==mini))
mini =
-0.1860
tmini =
-1.2850
6 ode45初试
format long
ts=[0,1];
y0=[1;0];
dydt=@(t,y)[y(2);-2*y(1)+3*y(2)+1]; %匿名函数写成的ode45所需dydt [tt,yy]=ode45(dydt,ts,y0);
y05=interp1(tt,yy(:,1),0.5,'spline')
%用一维插值求y(0.5)
y05 =
0.789580207901271
%符号运算法
syms t;
ys=dsolve('D2y-3*Dy+2*y=1','y(0)=1,Dy(0)=0','t') ;
ys_05=subs(ys,t,sym('0.5'))
ys_05 =
0.78958035647060552916850705213783
7 值空间基阵(第一天参考)
A=magic(8);
format short;
[R,jb]=rref(A);
B1=A(:,jb) %基向量
B2=orth(A) %单位正交基
B1 =
64 2 3
9 55 54
17 47 46
40 26 27
32 34 35
41 23 22
49 15 14
8 58 59
B2 =
-0.3536 0.5401 0.3536
-0.3536 -0.3858 -0.3536
-0.3536 -0.2315 -0.3536
-0.3536 0.0772 0.3536
-0.3536 -0.0772 0.3536
-0.3536 0.2315 -0.3536
-0.3536 0.3858 -0.3536
-0.3536 -0.5401 0.3536
B1_A=rref([B1,A])
B2_A=rref([B2,A])
B1_A =
1 0 0 1 0 0 1 1 0 0 1
0 1 0 0 1 0 3 4 -3 -4 7
0 0 1 0 0 1 -3 -4 4 5 -7
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
B2_A =
Columns 1 through 10
1.0000 0 0 -91.9239 -91.9239 -91.9239 -91.9239
-91.9239 -91.9239 -91.9239
0 1.0000 0 51.8459 -51.8459 -51.8459 51.8459
51.8459 -51.8459 -51.8459
0 0 1.0000 9.8995 -7.0711 -4.2426 1.4142
-1.4142 4.2426 7.0711
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
Column 11
-91.9239
51.8459
-9.8995
8 gallery矩阵
A=gallery(5)
[V,D]=eig(A);
diag(D)'
A =
-9 11 -21 63 -252
70 -69 141 -421 1684
-575 575 -1149 3451 -13801
3891 -3891 7782 -23345 93365
1024 -1024 2048 -6144 24572
ans =
-0.0405 + 0.0000i -0.0118 - 0.0383i -0.0118 + 0.0383i 0.0320 - 0.0228i 0.0320 + 0.0228i
%验算看分解精度
Areset=V*D/V
eAreset=norm(A-Areset,'fro')/norm(A,'fro')
Areset =
1.0e+04 *
-0.0009 - 0.0000i 0.0011 + 0.0000i -0.0021 - 0.0000i 0.0063 + 0.0000i -0.0252 - 0.0000i
0.0070 + 0.0000i -0.0069 - 0.0000i 0.0141 + 0.0000i -0.0421 - 0.0000i 0.1684 + 0.0000i
-0.0575 - 0.0000i 0.0575 + 0.0000i -0.1149 - 0.0000i 0.3451 + 0.0000i -1.3801 - 0.0000i
0.3891 + 0.0000i -0.3891 - 0.0000i 0.7782 + 0.0000i -2.3345 - 0.0000i 9.3365 + 0.0000i
0.1024 + 0.0000i -0.1024 - 0.0000i 0.2048 + 0.0000i -0.6144 - 0.0000i 2.4572 + 0.0000i
eAreset =
3.0594e-06
精度在1e-6左右,可以接受。
[V,D,s]=condeig(A)
V =
-0.0000 + 0.0000i -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i 0.0000 - 0.0000i
0.0206 + 0.0000i 0.0206 + 0.0001i 0.0206 - 0.0001i 0.0207 + 0.0001i 0.0207 - 0.0001i
-0.1398 + 0.0000i -0.1397 + 0.0001i -0.1397 - 0.0001i -0.1397 + 0.0000i -0.1397 - 0.0000i
0.9574 + 0.0000i 0.9574 + 0.0000i 0.9574 + 0.0000i 0.9574 + 0.0000i 0.9574 + 0.0000i
0.2519 + 0.0000i 0.2519 - 0.0000i 0.2519 + 0.0000i 0.2519 - 0.0000i
0.2519 + 0.0000i
D =
-0.0405 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i -0.0118 + 0.0383i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i -0.0118 - 0.0383i 0.0000 + 0.0000i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0320 + 0.0228i 0.0000 + 0.0000i
0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0000 + 0.0000i 0.0320 - 0.0228i
s =
1.0e+10 *
2.1969
2.1468
2.1468
2.0688
2.0688
其中,当s较大时,对应的特征值可信度较差
9 求解矩阵
A=magic(3);
b=ones(3,1)
x1=A\b
b =
1
1
1
x1 =
0.0667
0.0667
0.0667
[R,C]=rref([A,b]) %消去法解方程,同时判断解的唯一性
R =
1.0000 0 0 0.0667
0 1.0000 0 0.0667
0 0 1.0000 0.0667
C =
1 2 3
x3=inv(A)*b %左乘A逆解方程
x3 =
0.0667
0.0667
0.0667
三种方法所得解形式一致
10 四阶魔方阵
A=magic(4);
b=ones(4,1);
x1=A\b
x2=inv(A)*b
[x3,C]=rref([A,b])
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-17.
x1 =
0.0588
0.1176
-0.0588
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-17.
x2 =
0.0469
0.1875
-0.0625
-0.0156
x3 =
1.0000 0 0 1.0000 0.0588
0 1.0000 0 3.0000 0.1176
0 0 1.0000 -3.0000 -0.0588
0 0 0 0 0
C =
1 2 3
可以观察到,A不满秩,用逆矩阵的方法无法求得正确的解。
b可以由A的前三列线性表出,解不唯一。
验算:
A1=A*x1
A2=A*x2
A3=A(:,C)*x3(C,5)
A1 =
1
1
1
1
A2 =
0.7344
1.5469
1.1719
1.8594
A3 =
1
1
1
1
验算表明,逆矩阵法结果有误,理论上逆矩阵不存在。
xx=null(A) %齐次解
xx =
0.2236
0.6708
-0.6708
-0.2236
f=randn(1,6); %6个随机系数
xf=repmat(x1,1,6)+xx*f; %产生6个不同的特解
A*xf
ans =
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
11 求解矩阵
A=magic(4);
b=(1:4)';
x1=A\b
x2=inv(A)*b
[x3,C]=rref([A,b])
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-17.
x1 =
1.0e+15 *
-0.5629
-1.6888
1.6888
0.5629
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.306145e-17.
x2 =
1.0e+15 *
-0.5629
-1.6888
1.6888
0.5629
x3 =
1 0 0 1 0
0 1 0 3 0
0 0 1 -3 0
0 0 0 0 1
C =
1 2 3 5
由于该矩阵不满秩,b也不再A的值空间中,因而无法求得准确解。
12 实数解
t=-5:0.01:5;
y=@(t)(-0.5)+t-10*exp(-0.2*t).*abs(sin(sin(t)));
plot(t,y(t)) %利用匿名函数求y函数值
grid on;
[tt1,yy1]=ginput(1) %从图形获得近似解
tt1 =
2.7535
yy1 =
0.0365
[t,yt]=fzero(y,tt1)%寻找函数零点,接近的地方t =
2.7341
yt =
-4.4409e-16
13 求解二元函数方程组
%符号法
S=solve('sin(x-y)','cos(x+y)','x','y'); X=S.x
Y=S.y
X =
(3*pi)/4
-(3*pi)/4
pi/4
-pi/4
Y =
-pi/4
pi/4
-(3*pi)/4
(3*pi)/4
%图像交点法
x=-5:0.1:5;
y=x;
[X,Y]=meshgrid(x,y);
z1=sin(X-Y);
z2=cos(X+Y);
axis square;
grid on;
hold on;
contour(X,Y,z1,3) contour(X,Y,z2,3)
14 概率分布曲线N=100;
p=0.157;
P=binopdf(0:N,N,p); stem(0:N,P);
grid on;
15 正态分布随机数组
rng default
mu=4;sigma=2;
a=normrnd(mu,sigma,10000,1); subplot(211);
hist(a);
subplot(212);
histfit(a);
16 程序正确吗?
load prob_data416;%R是一个二维数组
Mx=max(max(R))
Me=mean(mean(R))
St=std(std(R))
Mx =
0.9997
Me =
0.5052
St =
0.0142
其中标准差是各列数据的标准差,与整个数组标准差有区别。
Mx1=max(R(:))
Me1=mean(R(:))
S1=std(R(:))
Mx1 =
0.9997
Me1 =
0.5052
S1 =
0.2919
其中标准差是整个数组标准差
17 求分式多项式,商多项式
N=conv([3 0 1 0],[1 0 0 0.5]);%用离散卷积表示多项式相乘
D=conv([1 2 -2],[5 2 0 1]);
[Q,r]=deconv(N,D)
Q =
0.6000 -1.4400
r =
0 0 21.8800 -5.3400 -5.5200 4.5800 -2.8800 yif=conv(D,Q)+r
yif =
3.0000 0 1.0000 1.5000 0 0.5000 0 err=norm(N-yif)
err =
7.1607e-15
数值计算存在误差,可以近似认为相等。
18 五阶拟合多项式
load prob_data418;
poly=polyfit(x,y,5);
yy=polyval(poly,x);
plot(x,y,'.b');
hold on;
plot(x,yy,'-r');
19 系统冲激响应
h=[0.05,0.24,0.40,0.24,0.15,-0.1,0.1]; rng default;
u=2*(randn(1,100)>0.5)-1;
con=conv(u,h);
lu=length(u);
lc=length(con);
u1=[u,nan*ones(1,lc-lu)];
subplot(211);
stem(0:lc-1,u1,'filled');
axis([0 lc -1 1]);
title('Input u');
subplot(212);
stem(0:lc-1,con,'filled');
axis([0 lc -1 1]);
title('Output y')
第五章
1 椭圆绘制
t=0:pi/100:2*pi;
a=4;b=2;
x=a*cos(t);y=b*sin(t);
plot(x,y,'.r','markersize',6) xlabel('x');ylabel('y');
axis equal;
2 心脏线
theta=0:pi/100:2*pi;
rou=1-cos(theta);
h=polar(theta,rou)
set(h,'linewidth',4,'color','r'); axis equal;
title('\rho =1-cos \theta')
直角坐标绘图
t=0:pi/100:2*pi;
x=2*cos(t)-cos(2*t);
y=2*sin(t)-sin(2*t);
plot(x,y,'r','linewidth',4)
axis equal;grid on;
title('\rho =1-cos \theta')
3 直方图
x=1:6;x=x';
y=[170,120,180,200,190,220;120,100,110,180,170,180;70,50,80,100,95,120] ';
bar(x,y,'stacked');
axis([0,7,0,600]);
colormap(cool);
legend('A','B','C','Location','NorthWest');
4 二阶线性系统的归一化
% exmp504.m 供第4道习题使用的程序
clc,clf,clear;
t=(0:0.05:18)';
N=length(t);
zeta=linspace(0.2,1.4,7);%改变产生向量的方式
L=length(zeta);
y=zeros(N,L);
hold on
for k=1:L
zk=zeta(k);
beta=sqrt(abs(1-zk^2));
if zk<1 %满足此条件,绘蓝色线
y=1/beta*exp(-zk*t).*sin(beta*t);
plot(t,y,'b')
if zk<0.4
text(2.2,0.63,'\zeta = 0.2')
end
elseif zk==1 %满足此条件,绘黑色线;问题在于可能无法取得 y=t.*exp(-t);
plot(t,y,'k','LineWidth',2)
zk
else %其余,绘红色线
y=(exp(-(zk-beta)*t)-exp(-(zk+beta)*t))/(2*beta);
plot(t,y,'r')
if zk>1.2
text(0.3,0.14,'\zeta = 1.4')
end
end
end
text(10,0.7,'\Delta\zeta=0.2')
axis([0,18,-0.4,0.8])
hold off
box on
grid on
zk =
1。