第四章matlab攻略.docx
- 文档编号:17431505
- 上传时间:2023-07-25
- 格式:DOCX
- 页数:50
- 大小:877.89KB
第四章matlab攻略.docx
《第四章matlab攻略.docx》由会员分享,可在线阅读,更多相关《第四章matlab攻略.docx(50页珍藏版)》请在冰点文库上搜索。
第四章matlab攻略
第4章数值计算
与符号计算相比,数值计算在科研和工程中的应用更为广泛。
MATLAB也正是凭借其卓越的数值计算能力而称雄世界。
随着科研领域、工程实践的数字化进程的深入,具有数字化本质的数值计算就显得愈益重要。
较之十年、二十年前,在计算机软硬件的支持下当今人们所能拥有的计算能力已经得到了巨大的提升。
这就自然地激发了人们从新的计算能力出发去学习、理解概念的欲望,鼓舞了人们用新计算能力试探解决实际问题的雄心。
鉴于当今高校本科教学偏重符号计算和便于手算简单示例的实际,也出于帮助读者克服对数值计算生疏感的考虑,本章在内容安排上仍从“微积分”开始。
这一方面与第2章符号计算相呼应,另方面通过“微积分”说明数值计算离散本质的微观和宏观影响。
为便于读者学习,本章内容的展开脉络基本上沿循高校数学教程,而内容深度力求控制在高校本科水平。
考虑到知识的跳跃和交叉,本书对重要概念、算式、指令进行了尽可能完整地说明。
.1数值微积分
.1.1近似数值极限及导数
【例4.1-1】设
,
,试用机器零阈值eps替代理论0计算极限
,
。
x=eps;
L1=(1-cos(2*x))/(x*sin(x)),
L2=sin(x)/x,
L1=
0
L2=
1
symst
f1=(1-cos(2*t))/(t*sin(t));
f2=sin(t)/t;
Ls1=limit(f1,t,0)
Ls2=limit(f2,t,0)
Ls1=
2
Ls2=
1
【例4.1-2】已知
,求该函数在区间
中的近似导函数。
d=pi/100;
t=0:
d:
2*pi;
x=sin(t);
dt=5*eps;
x_eps=sin(t+dt);
dxdt_eps=(x_eps-x)/dt;
plot(t,x,'LineWidth',5)
holdon
plot(t,dxdt_eps)
holdoff
legend('x(t)','dx/dt')
xlabel('t')
图4.1-1增量过小引起有效数字严重丢失后的毛刺曲线
x_d=sin(t+d);
dxdt_d=(x_d-x)/d;
plot(t,x,'LineWidth',5)
holdon
plot(t,dxdt_d)
holdoff
legend('x(t)','dx/dt')
xlabel('t')
图4.1-2增量适当所得导函数比较光滑
【例4.1-3】已知
,采用diff和gradient计算该函数在区间
中的近似导函数。
。
clf
d=pi/100;
t=0:
d:
2*pi;
x=sin(t);
dxdt_diff=diff(x)/d;
dxdt_grad=gradient(x)/d;
subplot(1,2,1)
plot(t,x,'b')
holdon
plot(t,dxdt_grad,'m','LineWidth',8)
plot(t(1:
end-1),dxdt_diff,'.k','MarkerSize',8)
axis([0,2*pi,-1.1,1.1])
title('[0,2\pi]')
legend('x(t)','dxdt_{grad}','dxdt_{diff}','Location','North')
xlabel('t'),boxoff
holdoff
subplot(1,2,2)
kk=(length(t)-10):
length(t);
holdon
plot(t(kk),dxdt_grad(kk),'om','MarkerSize',8)
plot(t(kk-1),dxdt_diff(kk-1),'.k','MarkerSize',8)
title('[end-10,end]')
legend('dxdt_{grad}','dxdt_{diff}','Location','SouthEast')
xlabel('t'),boxoff
holdoff
图4.1-3diff和gradient求数值近似导数的异同比较
.1.2数值求和与近似数值积分
【例4.1-4】求积分
,其中
。
clear
d=pi/8;
t=0:
d:
pi/2;
y=0.2+sin(t);
s=sum(y);
s_sa=d*s;
s_ta=d*trapz(y);
disp(['sum求得积分',blanks(3),'trapz求得积分'])
disp([s_sa,s_ta])
t2=[t,t(end)+d];
y2=[y,nan];
stairs(t2,y2,':
k')
holdon
plot(t,y,'r','LineWidth',3)
h=stem(t,y,'LineWidth',2);
set(h
(1),'MarkerSize',10)
axis([0,pi/2+d,0,1.5])
holdoff
shg
sum求得积分trapz求得积分
1.57621.3013
图4.1-4sum和trapz求积模式示意
.1.3计算精度可控的数值积分
【例4.1-5】求
。
symsx
Isym=vpa(int(exp(-x^2),x,0,1))
Isym=
.74682413281242702539946743613185
formatlong
d=0.001;x=0:
d:
1;
Itrapz=d*trapz(exp(-x.*x))
Itrapz=
0.74682407149919
fx='exp(-x.^2)';
Ic=quad(fx,0,1,1e-8)
Ic=
0.74682413285445
【例4.1-6】求
。
symsxy
s=vpa(int(int(x^y,x,0,1),y,1,2))
s=
.40546510810816438197801311546432
formatlong
s_n=dblquad(@(x,y)x.^y,0,1,1,2)
s_n=
0.40546626724351
.1.4函数极值的数值求解
【例4.1-7】已知
,在
区间,求函数的极小值。
(1)
symsx
y=(x+pi)*exp(abs(sin(x+pi)));
yd=diff(y,x);
xs0=solve(yd)
yd_xs0=vpa(subs(yd,x,xs0),6)
y_xs0=vpa(subs(y,x,xs0),6)
y_m_pi=vpa(subs(y,x,-pi/2),6)
y_p_pi=vpa(subs(y,x,pi/2),6)
xs0=
-1.0676598444985783372948670854801
yd_xs0=
.1e-4
y_xs0=
4.98043
y_m_pi=
4.26987
y_p_pi=
12.8096
(2)
x1=-pi/2;x2=pi/2;
yx=@(x)(x+pi)*exp(abs(sin(x+pi)));
[xn0,fval,exitflag,output]=fminbnd(yx,x1,x2)
xn0=
-1.2999e-005
fval=
3.1416
exitflag=
1
output=
iterations:
21
funcCount:
24
algorithm:
'goldensectionsearch,parabolicinterpolation'
message:
[1x112char]
(3)
xx=-pi/2:
pi/200:
pi/2;
yxx=(xx+pi).*exp(abs(sin(xx+pi)));
plot(xx,yxx)
xlabel('x'),gridon
图4.1-5在[-pi/2,pi/2]区间中的函数曲线
[xx,yy]=ginput
(1)
xx=
1.5054e-008
yy=
3.1416
图4.1-6函数极值点附近的局部放大和交互式取值
【例4.1-8】求
的极小值点。
它即是著名的Rosenbrock's"Banana"测试函数,它的理论极小值是
。
(1)
ff=@(x)(100*(x
(2)-x
(1).^2)^2+(1-x
(1))^2);
(2)
x0=[-5,-2,2,5;-5,-2,2,5];
[sx,sfval,sexit,soutput]=fminsearch(ff,x0)
sx=
1.0000-0.68970.41518.0886
1.0000-1.91684.96437.8004
sfval=
2.4112e-010
sexit=
1
soutput=
iterations:
384
funcCount:
615
algorithm:
'Nelder-Meadsimplexdirectsearch'
message:
[1x196char]
(3)检查目标函数值
formatshorte
disp([ff(sx(:
1)),ff(sx(:
2)),ff(sx(:
3)),ff(sx(:
4))])
2.4112e-0105.7525e+0022.2967e+0033.3211e+005
.1.5常微分方程的数值解
【例4.1-9】求微分方程
,
,在初始条件
情况下的解,并图示。
(见图4.1-7和4.1-8)
(1)
(2)
[DyDt.m]
functionydot=DyDt(t,y)
mu=2;
ydot=[y
(2);mu*(1-y
(1)^2)*y
(2)-y
(1)];
(3)解算微分方程
tspan=[0,30];
y0=[1;0];
[tt,yy]=ode45(@DyDt,tspan,y0);
plot(tt,yy(:
1))
xlabel('t'),title('x(t)')
图4.1-7微分方程解
(4)
plot(yy(:
1),yy(:
2))
xlabel('位移'),ylabel('速度')
图4.1-8平面相轨迹
.2矩阵和代数方程
.2.1矩阵运算和特征参数
10一矩阵运算
【例4.2-1】已知矩阵
,
,采用三种不同的编程求这两个矩阵的乘积
。
(1)
clear
rand('state',12)
A=rand(2,4);B=rand(4,3);
%------------------
C1=zeros(size(A,1),size(B,2));
forii=1:
size(A,1)
forjj=1:
size(B,2)
fork=1:
size(A,2)
C1(ii,jj)=C1(ii,jj)+A(ii,k)*B(k,jj);
end
end
end
C1
C1=
1.05590.98591.2761
1.76541.84831.8811
(2)
%------------------
C2=zeros(size(A,1),size(B,2));
forjj=1:
size(B,2)
fork=1:
size(B,1)
C2(:
jj)=C2(:
jj)+A(:
k)*B(k,jj);
end
end
C2
C2=
1.05590.98591.2761
1.76541.84831.8811
(3)
C3=A*B,
C3=
1.05590.98591.2761
1.76541.84831.8811
(3)
C3_C3=norm(C3-C1,'fro')
C3_C2=norm(C3-C2,'fro')
C3_C3=
0
C3_C2=
0
【例4.2-2】观察矩阵的转置操作和数组转置操作的差别。
formatrat
A=magic
(2)+j*pascal
(2)
A=
1+1i3+1i
4+1i2+2i
%
A1=A'
A2=A.'
A1=
1-1i4-1i
3-1i2-2i
A2=
1+1i4+1i
3+1i2+2i
B1=A*A'
B2=A.*A'
C1=A*A.'
C2=A.*A.'
B1=
1213-1i
13+1i25
B2=
213+1i
13-1i8
C1=
8+8i7+13i
7+13i15+16i
C2=
0+2i11+7i
11+7i0+8i
10二矩阵的标量特征参数
【例4.2-3】矩阵标量特征参数计算示例。
A=reshape(1:
9,3,3)
r=rank(A)
d3=det(A)
d2=det(A(1:
2,1:
2))
t=trace(A)
A=
147
258
369
r=
2
d3=
0
d2=
-3
t=
15
【例4.2-4】矩阵标量特征参数的性质。
formatshortg
rand('state',0)
A=rand(3,3);
B=rand(3,3);
C=rand(3,4);
D=rand(4,3);
tAB=trace(A*B)
tBA=trace(B*A)
tCD=trace(C*D)
tDC=trace(D*C)
tAB=
3.6697
tBA=
3.6697
tCD=
2.1544
tDC=
2.1544
d_A_B=det(A)*det(B)
dAB=det(A*B)
dBA=det(B*A)
d_A_B=
0.0846
dAB=
0.0846
dBA=
0.0846
dCD=det(C*D)
dDC=det(D*C)
dCD=
-0.012362
dDC=
-2.1145e-018
.2.2矩阵的变换和特征值分解
【例4.2-5】行阶梯阵简化指令rref计算结果的含义。
(1)
A=magic(4)
[R,ci]=rref(A)
A=
162313
511108
97612
414151
R=
1001
0103
001-3
0000
ci=
123
(2)
r_A=length(ci)
r_A=
3
(3)
aa=A(:
1:
3)*R(1:
3,4)
err=norm(A(:
4)-aa)
aa=
13
8
12
1
err=
0
【例4.2-6】矩阵零空间及其含义。
A=reshape(1:
15,5,3);
X=null(A)
S=A*X
n=size(A,2);
l=size(X,2);
n-l==rank(A)
X=
-0.4082
0.8165
-0.4082
S=
1.0e-014*
0.2109
0.2220
0.3220
0.3331
0.4330
ans=
1
【例4.2-7】简单实阵的特征值分解。
(1)
A=[1,-3;2,2/3]
[V,D]=eig(A)
A=
1.0000-3.0000
2.00000.6667
V=
0.77460.7746
0.0430-0.6310i0.0430+0.6310i
D=
0.8333+2.4438i0
00.8333-2.4438i
(2)
[VR,DR]=cdf2rdf(V,D)
VR=
0.77460
0.0430-0.6310
DR=
0.83332.4438
-2.44380.8333
(3)
A1=V*D/V
A1_1=real(A1)
A2=VR*DR/VR
err1=norm(A-A1,'fro')
err2=norm(A-A2,'fro')
A1=
1.0000-0.0000i-3.0000
2.0000+0.0000i0.6667
A1_1=
1.0000-3.0000
2.00000.6667
A2=
1.0000-3.0000
2.00000.6667
err1=
4.4409e-016
err2=
4.4409e-016
.2.3线性方程的解
10一线性方程解的一般结论
10二除法运算解方程
【例4.2-8】求方程
的解。
(1)
A=reshape(1:
12,4,3);
b=(13:
16)';
(2)
ra=rank(A)
rab=rank([A,b])
ra=
2
rab=
2
(3)
xs=A\b;
xg=null(A);
c=rand
(1);
ba=A*(xs+c*xg)
norm(ba-b)
Warning:
Rankdeficient,rank=2,tol=1.8757e-014.
ba=
13.0000
14.0000
15.0000
16.0000
ans=
1.3874e-014
10三矩阵逆
【例4.2-9】“逆阵”法和“左除”法解恰定方程的性能对比
(1)
randn('state',0);
A=gallery('randsvd',300,2e13,2);
x=ones(300,1);
b=A*x;
cond(A)
ans=
2.0069e+013
(2)
tic
xi=inv(A)*b;
ti=toc
eri=norm(x-xi)
rei=norm(A*xi-b)/norm(b)
ti=
0.0341
eri=
0.0839
rei=
0.0048
(3)
tic;
xd=A\b;
td=toc
erd=norm(x-xd)
red=norm(A*xd-b)/norm(b)
td=
0.0172
erd=
0.0127
red=
9.4744e-015
.2.4一般代数方程的解
【例4.2-10】求
的零点。
(1)
S=solve('sin(t)^2*exp(-0.1*t)-0.5*abs(t)','t')
S=
0.
(2)
y_C=inline('sin(t).^2.*exp(-0.1*t)-0.5*abs(t)','t');
t=-10:
0.01:
10;
Y=y_C(t);
clf,
plot(t,Y,'r');
holdon
plot(t,zeros(size(t)),'k');
xlabel('t');ylabel('y(t)')
holdoff
图4.2-1函数零点分布观察图
zoomon
[tt,yy]=ginput(5);zoomoff
图4.2-2局部放大和利用鼠标取值图
tt
tt=
-2.0039
-0.5184
-0.0042
0.6052
1.6717
[t4,y4]=fzero(y_C,0.1)
t4=
0.5993
y4=
1.1102e-016
.3概率分布和统计分析
.3.1概率函数、分布函数、逆分布函数和随机数的发生
10一二项分布(Binomialdistribution)
【例4.3-1】画出N=100,p=0.5情况下的二项分布概率特性曲线。
N=100;p=0.5;
k=0:
N;
pdf=binopdf(k,N,p);
cdf=binocdf(k,N,p);
h=plotyy(k,pdf,k,cdf);
set(get(h
(1),'Children'),'Color','b','Marker','.','MarkerSize',13)
set(get(h
(1),'Ylabel'),'String','pdf')
set(h
(2),'Ycolor',[1,0,0])
set(get(h
(2),'Children'),'Color','r','Marker','+','MarkerSize',4)
set(get(h
(2),'Ylabel'),'String','cdf')
xlabel('k')
gridon
图4.3-1二项分布B(100,0.5)的概率和累计概率曲线
10二正态分布(Normaldistribution)
【例4.3-2】正态分布标准差的几何表示。
mu=3;sigma=0.5;
x=mu+sigma*[-3:
-1,1:
3];
yf=normcdf(x,mu,sigma);
P=[yf(4)-yf(3),yf(5)-yf
(2),yf(6)-yf
(1)];
xd=1:
0.1:
5;
yd=normpdf(xd,mu,sigma);
clf
fork=1:
3
%-------------------------------
xx=x(4-k):
sigma/10:
x(3+k);
yy=normpdf(xx,mu,sigma);
%--------------------------------
subplot(3,1,k),plot(xd,yd,'b');
holdon
fill([x(4-k),xx,x(3+k)],[0,yy,0],'g');
holdoff
ifk<2
text(3.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 第四 matlab 攻略