ch4数值计算.docx

ch4数值计算.docx

《ch4数值计算.docx》由会员分享，可在线阅读，更多相关《ch4数值计算.docx（48页珍藏版）》请在冰点文库上搜索。

ch4数值计算

第4章数值计算

.1数值微积分

.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）

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]）

legend（'g_1（x）=1/x^{0.2}','g_2（x）=x^5','Location','North'）

title（'x位于[0,1]间的g_1（x）曲线与g_2（x）曲线所夹的区域'）

图4.1-5待求面积的区域

（2）采用标量型匿名函数计算积分

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%<12>

Q1=

1.250000027856048

Q2=

0.166********6667

S12=

.0833********

（3）

G=@（x）x.^[-0.2;5];%<13>

Q=integral（G,0,1,'ArrayValued',true）%<14>

S=[1,-1]*Q%<15>

Q=

1.250000027856048

0.166********6667

S=

1.083333361189381

（4）

symst%

Gsym=vpa（int（t.^[-0.2;5],0,1））;%

Ssym=Gsym

（1）-Gsym

（2）%<17>

Ssym=

1.0833333333333333333333333333333

【例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,'Waypoints',Path）

sf=

0.0000+12.5664i

.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=

0.050838341410271656880659496266968

yd_xs0=

2.29589*10^（-41）

y_xs0=

-0.00126332

（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矩阵的标量特征参数

【例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=

3.7479

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=

-0.0852

dBA=

-0.0852

dCD=det（C*D）

dDC=det（D*C）%

dCD=

-0.0557

dDC=

1.5085e-017

.2.2矩阵的变换和特征值分解

【例4.2-3】

（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-4】

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.2665

0.2665

0.3553

0.4441

0.6217

ans=

1

【例4.2-5】

（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-3.0000-0.0000i

2.00000.6667

A1_1=

1.0000-3.0000

2.00000.6667

A2=

1.0000-3.0000

2.00000.6667

err1=

4.6613e-016

err2=

4.4409e-016

.2.3线性方程的解

101线性方程解的一般结论

102除法运算解方程

【例4.2-6】

（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

103矩阵逆

【例4.2-7】

（1）

rngdefault

A=gallery（'randsvd',300,2e13,2）;%

x=ones（300,1）;%

b=A*x;%

cond（A）%

ans=

2.0022e+013

（2）

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

（3）

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一般代数方程的解

【例4.2-8】

（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=

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】

图4.3-3概率分布交互界面

.3.2全局随机流、随机数组和统计分析

101全局随机流的操控

表4.3-1产生全局随机流的发生器

generator

算例

可取字符串

含义

'twister'

默认的随机数发生器

例

'v5uniform'

MATLAB5.0版均布随机数发生器

'v5normal'

MATLAB5.0版正态随机数发生

例4.3-5,例4.3-6

'v4'

MATLAB4.0版随机发生器

102三个基本随机数组创建指令

【例4.3-4】

（1）

rngdefault%<1>

GRS=rand（1,25）%<2>

GRS=

Columns1through5

0.81470.90580.12700.91340.6324

Columns6through10

0.09750.27850.54690.95750.9649

Columns11through15

0.15760.97060.95720.48540.8003

Columns16through20

0.14190.42180.91570.79220.9595

Columns21through25

0.65570.03570.84910.93400.6787

（2）

rngdefault%<3>

r1=randn（1,5）%<4>

r2=rand（1,5）%<5>

r3=randi（[-3,2],1,5）%<6>

r4