ch4数值计算Word文档格式.docx
《ch4数值计算Word文档格式.docx》由会员分享，可在线阅读，更多相关《ch4数值计算Word文档格式.docx（46页珍藏版）》请在冰豆网上搜索。
holdon
plot（t,dxdt_d）
）
图4.1-2
【例4.1-3】
clf
dxdt_diff=diff（x）/d;
dxdt_grad=gradient（x）/d;
%
subplot（1,2,1）
b'
plot（t,dxdt_grad,'
m'
8）
plot（t（1:
end-1）,dxdt_diff,'
.k'
MarkerSize'
axis（[0,2*pi,-1.1,1.1]）
title（'
[0,2\pi]'
dxdt_{grad}'
dxdt_{diff}'
Location'
North'
）,boxoff
subplot（1,2,2）
kk=（length（t）-10）:
length（t）;
plot（t（kk）,dxdt_grad（kk）,'
om'
plot（t（kk-1）,dxdt_diff（kk-1）,'
[end-10,end]'
SouthEast'
holdoff
图4.1-3
.1.2数值求和与近似数值积分
【例4.1-4】
clear
d=pi/8;
%
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];
hs=stairs（t2,y2,'
:
k'
3）;
%<
12>
ht=plot（t,y,'
r'
%<
14>
stem（t,y）%
legend（[hs,ht],'
sum'
trapz'
location'
best'
）<
16>
axis（[0,pi/2+d,0,1.5]）%
shg
sum求得积分trapz求得积分
1.57621.3013
图4.1-4
.1.3计算精度可控的数值积分
表4.2-1积分指令integral的选项名称、取值、及默认值
选项名
Name
允许值
Value
默认值
DefaultValue
含义
'
AbsTol'
正实数
RelTol'
ArrayValued'
false
true
Waypoints'
实数单调序列数组
复数序列数组
【例4.1-5】在
之间，求曲线
与
所夹区域的面积（参见图4.1-5）。
x=linspace（0.01,1.2,50）;
g1=x.^（-0.2）;
g2=x.^5;
plot（x,g1,'
-r.'
x,g2,'
-b*'
axis（[0,1.2,0,3]）
g_1（x）=1/x^{0.2}'
g_2（x）=x^5'
x位于[0,1]间的g_1（x）曲线与g_2（x）曲线所夹的区域'
）
图4.1-5
formatlong
G1=@（x）x.^-0.2;
8>
Q1=integral（G1,0,1）%<
9>
G2=@（x）x.^5;
%<
10>
Q2=integral（G2,0,1）%<
11>
S12=Q1-Q2%<
Q1=
1.250000027856048
Q2=
0.166********6667
S12=
.0833********
3）
G=@（x）x.^[-0.2;
5];
%<
13>
Q=integral（G,0,1,'
true）%<
S=[1,-1]*Q%<
15>
Q=
S=
1.083333361189381
4）
symst%
Gsym=vpa（int（t.^[-0.2;
5],0,1））;
Ssym=Gsym
（1）-Gsym
（2）%
Ssym=
1.0833333333333333333333333333333
【例4.1-6】验证
图4.1-6
F=@（z）（2*z-1）./（z.*（z-1））;
1>
Path=[2+1i,-1+1i,-1-1i,2-1i,2+1i];
2>
%
sf=integral（F,2+1i,2+1i,'
Path）
sf=
0.000000000000000+12.566370614359174i
.1.4函数极值的数值求解
【例4.1-7】
x1=-50;
x2=5;
yx=@（x）（sin（x）.^2.*exp（-0.1*x）-0.5*sin（x）.*（x+0.1））;
%<
%
[xc0,fc0,exitflag]=fminbnd（yx,x1,x2）%<
3>
%
xc0=
-8.4867
fc0=
-1.8621
exitflag=
xx=-50:
pi/200:
5;
ezplot（yx,[-50,5]）%<
5>
x'
）,gridon
xx=[-23,-20,-18];
fc=fc0;
xc=xc0;
fork=1:
[xw,fw]=fminbnd（yx,xx（k）,xx（k+1））;
%<
iffw<
fc0
xc=xw;
fc=fw;
end
end
fprintf（'
函数最小值%6.5f发生在x=%6.5f处'
fc,xc）
函数最小值-3.34765发生在x=-19.60721处
【例4.1-8】求
的极小值点。
ff=@（x）（100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2）;
formatshortg
x0=[-5,-2,2,5;
-5,-2,2,5];
[sx,sfval,sexit,soutput]=fminsearch（ff,x0）
%
sx=
0.99998-0.689710.415078.0886
0.99997-1.91684.96437.8004
sfval=
2.4112e-10
sexit=
1
soutput=
iterations:
384
funcCount:
615
algorithm:
'
Nelder-Meadsimplexdirectsearch'
message:
[1x194char]
formatshorte
disp（[ff（sx（:
1））,ff（sx（:
2））,ff（sx（:
3））,ff（sx（:
4））]）
2.4112e-105.7525e+022.2967e+033.3211e+05
.1.5常微分方程的数值解
【例4.1-9】求微分方程
，在初始条件
情况下的解，并图示。
functionydot=DyDt（t,y）
mu=2;
ydot=[y
（2）;
mu*（1-y
（1）^2）*y
（2）-y
（1）];
%
tspan=[0,30];
y0=[1;
0];
[tt,yy]=ode45（@DyDt,tspan,y0）;
%<
plot（tt,yy（:
1））
）,title（'
plot（yy（:
1）,yy（:
2））%
位移'
）,ylabel（'
速度'
图4.1-7
.2矩阵和代数方程
.2.1矩阵的标量特征参数
表4.2-2
术语
数学含义
MATLAB指令
秩
Rank
rank（A）
迹
Trace
trace（A）
行列式Determinant
det（A）
【例4.2-1】
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-2】
formatshort%
rngdefault%
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.7479
tBA=
tCD=
3.3399
tDC=
3.3399
d_A_B=det（A）*det（B）
dAB=det（A*B）
dBA=det（B*A）
d_A_B=
-0.0852
dAB=
dBA=
-0.0852
dCD=det（C*D）
dDC=det（D*C）%
dCD=
-0.0557
dDC=
1.5085e-017
.2.2矩阵的变换和特征值分解
【例4.2-3】
A=magic（4）%
[R,ci]=rref（A）%
162313
511108
97612
414151
R=
1001
0103
001-3
0000
ci=
123
2）
r_A=length（ci）
r_A=
3
aa=A（:
1:
3）*R（1:
3,4）%
err=norm（A（:
4）-aa）%
aa=
13
8
12
err=
0
【例4.2-4】
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
1.0e-014*
0.2665
0.3553
0.4441
0.6217
ans=
【例4.2-5】
A=[1,-3;
2,2/3]
[V,D]=eig（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
[VR,DR]=cdf2rdf（V,D）
VR=
0.77460
0.0430-0.6310
DR=
0.83332.4438
-2.44380.8333
A1=V*D/V%
A1_1=real（A1）%
A2=VR*DR/VR
err1=norm（A-A1,'
fro'
err2=norm（A-A2,'
A1=
1.0000-3.0000-0.0000i
2.00000.6667
A1_1=
A2=
err1=
4.6613e-016
err2=
4.4409e-016
.2.3线性方程的解
101线性方程解的一般结论
102除法运算解方程
【例4.2-6】求方程
的解。
12,4,3）;
b=（13:
16）'
;
%
ra=rank（A）%
rab=rank（[A,b]）%
ra=
2
rab=
2
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
1.1374e-014
103矩阵逆
【例4.2-7】
rngdefault
A=gallery（'
randsvd'
300,2e13,2）;
x=ones（300,1）;
b=A*x;
cond（A）%
2.0022e+013
tic%
xi=inv（A）*b;
ti=toc%
eri=norm（x-xi）%
rei=norm（A*xi-b）/norm（b）%
ti=
0.2034
eri=
0.0812
rei=
0.0047
tic;
xd=A\b;
td=toc
erd=norm（x-xd）
red=norm（A*xd-b）/norm（b）
td=
0.0125
erd=
0.0274
red=
7.5585e-015
.2.4一般代数方程的解
解题步骤大致如下：
（1）
（2）
【例4.2-8】
ft=sin（t）^2*exp（-0.1*t）-0.5*abs（t）;
S=solve（ft,t）%<
ftS=subs（ft,t,S）%
ftS=
0
y_C=inline（'
sin（t）.^2.*exp（-0.1*t）-0.5*abs（t）'
）;
%
t=-10:
0.01:
10;
Y=y_C（t）;
clf,
plot（t,Y,'
plot（t,zeros（size（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）%<
18>
t4=
0.5993
y4=
1.1102e-016
.3概率分布和统计分析
.3.1概率函数、分布函数、逆分布函数和随机数的发生
101二项分布（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'
Marker'
.'
13）
Ylabel'
String'
pdf'
set（h
（2）,'
Ycolor'
[1,0,0]）
set（get（h
（2）,'
+'
4）
cdf'
%<
）%
gridon
图4.3-1
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:
yd=normpdf（xd,mu,sigma）;
3
%-------------------------------
xx=x（4-k）:
sigma/10:
x（3+k）;
yy=normpdf（xx,mu,sigma）;
%--------------------------------------------------
subplot（3,1,k）,plot（xd,yd,'
holdon
fill（[x（4-k）,xx,x（3+k）],[0,yy,0],'
g'
holdoff
ifk<
text（3.8,0.6,'
[{\mu}-{\sigma},{\mu}+{\sigma}]'
else
kk=int2str（k）;
text（3.8,0.6,['
[{\mu}-'
kk,'
{\sigma},{\mu}+'
{\sigma}]'
text（2.8,0.3,num2str（P（k）））;
shg
end
shg
图4.3-2
103各种概率分布的交互式观察界面
【例4.3-3】概率分布交互界面使用方法介绍。
图4.3-3
.3.2全局随机流、随机数组和统计分析
101全局随机流的操控
表4.3-1
generator
算例
可取字符串
含义
twister'
v5uniform'
v5normal'
v4'
102三个基本随机数组创建指令
【例4.3-4】
rngdefault
升级会员 