ch4数值计算a.docx
- 文档编号:30645002
- 上传时间:2023-08-18
- 格式:DOCX
- 页数:44
- 大小:268.04KB
ch4数值计算a.docx
《ch4数值计算a.docx》由会员分享,可在线阅读,更多相关《ch4数值计算a.docx(44页珍藏版)》请在冰豆网上搜索。
ch4数值计算a
第4章数值计算
.1.1近似数值极限及导数
【例4.1-1】
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】
(1)
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增量过小引起有效数字严重丢失后的毛刺曲线
(2)
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】
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);%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;%<6>
s_ta=d*trapz(y);%<7>
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】
(1)
symsx
Isym=vpa(int(exp(-x^2),x,0,1))
Isym=
0.74682413281242702539946743613185
(2)
formatlong%
d=0.001;x=0:
d:
1;
Itrapz=d*trapz(exp(-x.*x))
Itrapz=
0.746824071499185
(3)
fx='exp(-x.^2)';%
Ic=quad(fx,0,1,1e-8)%<7>
Ic=
0.746824132854452
【例4.1-6】
(1)
symsxy
s=vpa(int(int(x^y,x,0,1),y,1,2))%<2>
Warning:
Explicitintegralcouldnotbefound.
s=
0.40546510810816438197801311546435
(2)
formatlong
s_n=dblquad(@(x,y)x.^y,0,1,1,2)%<4>
s_n=
0.405466267243508
.1.4函数极值的数值求解
【例4.1-7】
(1)
symsx
y=sin(x)^2*exp(-0.1*x)-0.5*sin(x)*(x+0.1);
yd=diff(y,x);%
xs0=solve(yd,x)%<4>
yd_xs0=vpa(subs(yd,x,xs0),6)%
y_xs0=vpa(subs(y,x,xs0),6)%
xs0=
matrix([[0.050838341410271656880659496266968]])
yd_xs0=
2.2958874039497802890014385492622*10^(-41)
y_xs0=
-0.001263317776974196724544154344118
(2)
x1=-10;x2=10;%
yx=@(x)(sin(x)^2*exp(-0.1*x)-0.5*sin(x)*(x+0.1));
%
[xn0,fval,exitflag,output]=fminbnd(yx,x1,x2)%<9>
%
xn0=
2.514797840754235
fval=
-0.499312445280039
exitflag=
1
output=
iterations:
13
funcCount:
14
algorithm:
'goldensectionsearch,parabolicinterpolation'
message:
[1x112char]
(3)
xx=-10:
pi/200:
10;%
yxx=subs(y,x,xx);
plot(xx,yxx)
xlabel('x'),gridon
图4.1-5在[-10,10]区间中的函数曲线
(4)
x11=6;x2=10;%
yx=@(x)(sin(x)^2*exp(-0.1*x)-0.5*sin(x)*(x+0.1));
%
[xn00,fval,exitflag,output]=fminbnd(yx,x11,x2)%<16>
%
xn00=
.023*********
fval=
-3.568014059128578
exitflag=
1
output=
iterations:
9
funcCount:
10
algorithm:
'goldensectionsearch,parabolicinterpolation'
message:
[1x112char]
【例4.1-8】
(1)
ff=@(x)(100*(x
(2)-x
(1).^2)^2+(1-x
(1))^2);
(2)
formatshortg
x0=[-5,-2,2,5;-5,-2,2,5];%提供4个搜索起点
[sx,sfval,sexit,soutput]=fminsearch(ff,x0)
%sx给出一组使优化函数值非减的局部极小点
sx=
0.99998-0.689710.415078.0886
0.99997-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】
(1)
(3)
tspan=[0,30];%
y0=[1;0];%
[tt,yy]=ode45(@DyDt,tspan,y0);%<3>
plot(tt,yy(:
1))
xlabel('t'),title('x(t)')
图4.1-6微分方程解
(4)
plot(yy(:
1),yy(:
2))%
xlabel('位移'),ylabel('速度')
图4.1-7平面相轨迹
.2矩阵和代数方程
.2.1矩阵运算和特征参数
101矩阵运算
【例4.2-1】
(1)
clear
rand('twister',12)
A=rand(2,4);B=rand(4,3);
C1=zeros(size(A,1),size(B,2));%<4>
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%<13>
C1
C1=
0.73370.83960.3689
1.06241.17341.3787
(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=
0.73370.83960.3689
1.06241.17341.3787
(3)
C3=A*B,%
C3=
0.73370.83960.3689
1.06241.17341.3787
(4)
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
102矩阵的标量特征参数
【例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】
formatshort%
rand('twister',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=
2.6030
tBA=
2.6030
tCD=
4.1191
tDC=
4.1191
d_A_B=det(A)*det(B)
dAB=det(A*B)
dBA=det(B*A)
d_A_B=
0.0094
dAB=
0.0094
dBA=
0.0094
dCD=det(C*D)
dDC=det(D*C)%
dCD=
0.0424
dDC=
-2.6800e-018
.2.2矩阵的变换和特征值分解
【例4.2-5】
(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
-0.1776
-0.2665
-0.3553
-0.5329
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=
6.7532e-016
err2=
4.4409e-016
.2.3线性方程的解
101线性方程解的一般结论
【例4.2-8】
(1)
A=reshape(1:
12,4,3);%
b=(13:
16)';%
(2)
ra=rank(A)%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.1374e-014
102矩阵逆
【例4.2-9】
(1)
randn('state',0);
A=gallery('randsvd',300,2e13,2);%
x=ones(300,1);%
b=A*x;%
cond(A)%
ans=
1.9978e+013
(2)
tic%
xi=inv(A)*b;%
ti=toc%
eri=norm(x-xi)%
rei=norm(A*xi-b)/norm(b)%
ti=
0.0294
eri=
0.1003
rei=
0.0053
(3)
tic;
xd=A\b;%
td=toc
erd=norm(x-xd)
red=norm(A*xd-b)/norm(b)
td=
0.0134
erd=
0.0939
red=
8.4835e-015
.2.4一般代数方程的解
【例4.2-10】
(1)
symst
ft=sin(t)^2*exp(-0.1*t)-0.5*abs(t);
S=solve(ft,t)%<3>
ftS=subs(ft,t,S)%
S=
matrix([[0]])
ftS=
0
(2)
(A)
y_C=inline('sin(t).^2.*exp(-0.1*t)-0.5*abs(t)','t');
%
(B)
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函数零点分布观察图
(C)
zoomon%
[tt,yy]=ginput(5);zoomoff%
图4.2-2局部放大和利用鼠标取值图
tt%
tt=
-2.0039
-0.5184
-0.0042
0.6052
1.6717
(D)
[t4,y4]=fzero(y_C,0.1)%<17>
t4=
0.5993
y4=
1.1102e-016
.3概率分布和统计分析
101二项分布(Binomialdistribution)
【例4.3-1】
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)的概率和累计概率曲线
102正态分布(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.8,0.6,'[{\mu}-{\sigma},{\mu}+{\sigma}]')
else
kk=int2str(k);
text(3.8,0.6,['[{\mu}-',kk,'{\sigma},{\mu}+',kk,'{\sigma}]'])
end
text(2.8,0.3,num2str(P(k)));shg
end
xlabel('x');shg
图4.3-2均值两侧一、二、三倍标准差之间的概率
103各种概率分布的交互式观察界面
【例4.3-3
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- ch4 数值 计算