ch4数值计算.docx

ch4数值计算.docx

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

ch4数值计算

第4章数值计算

与符号计算相比，数值计算在科研和工程中的应用更为广泛。

MATLAB也正是凭借其卓越的数值计算能力而称雄世界。

随着科研领域、工程实践的数字化进程的深入，具有数字化本质的数值计算就显得愈益重要。

较之十年、二十年前，在计算机软硬件的支持下当今人们所能拥有的计算能力已经得到了巨大的提升。

这就自然地激发了人们从新的计算能力出发去学习、理解概念的欲望，鼓舞了人们用新计算能力试探解决实际问题的雄心。

鉴于当今高校本科教学偏重符号计算和便于手算简单示例的实际，也出于帮助读者克服对数值计算生疏感的考虑，本章在内容安排上仍从“微积分”开始。

这一方面与第2章符号计算相呼应，另方面通过“微积分”说明数值计算离散本质的微观和宏观影响。

为便于读者学习，本章内容的展开脉络基本上沿循高校数学教程，而内容深度力求控制在高校本科水平。

考虑到知识的跳跃和交叉，本书对重要概念、算式、指令进行了尽可能完整地说明。

.1数值微积分

.1.1近似数值极限及导数

〖说明〗

【例4.1-1】设

，

，试用最小正数realmin替代理论0计算极限

，

。

%

x=realmin;

L1=（1-cos（2*x））/（x*sin（x））,%

L2=sin（x）/x,%

L1=

NaN

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）;%

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-3

〖说明〗

.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];%

hs=stairs（t2,y2,':

k','LineWidth',3）;%<12>

holdon

ht=plot（t,y,'r','LineWidth',3）;%<14>

stem（t,y）%

legend（[hs,ht],'sum','trapz','location','best'）<16>

axis（[0,pi/2+d,0,1.5]）%

holdoff

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

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

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

sf=

0.000000000000000+12.566370614359174i

〖说明〗

.1.4函数极值的数值求解

〖说明〗

【例4.1-7】

1）

x1=-50;x2=5;%

yx=@（x）（sin（x）.^2.*exp（-0.1*x）-0.5*sin（x）.*（x+0.1））;%<2>

%

[xc0,fc0,exitflag]=fminbnd（yx,x1,x2）%<3>

%

xc0=

-8.4867

fc0=

-1.8621

exitflag=

1

2）

xx=-50:

pi/200:

5;%

ezplot（yx,[-50,5]）%<5>

xlabel（'x'）,gridon

图4.1-5

3）

xx=[-23,-20,-18];%

fc=fc0;xc=xc0;%

fork=1:

2

[xw,fw]=fminbnd（yx,xx（k）,xx（k+1））;%<16>

iffw

xc=xw;

fc=fw;

end

end

fprintf（'函数最小值%6.5f发生在x=%6.5f处',fc,xc）

函数最小值-3.34765发生在x=-19.60721处

〖说明〗

【例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];%

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

3）

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】求微分方程

，

，在初始条件

，

情况下的解，并图示。

1）

2）

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）;%<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-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=

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

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

解题步骤大致如下：

（1）

（2）

〖说明〗

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

%

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）%<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）;%<5>

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

%<10>

xlabel（'k'）%

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:

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】概率分布交互界面使用方法介绍。

1）

2）

3）

4）

图4.3-3

〖说明〗

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

101全局随机流的操控

〖说明〗

表4.3-1

generator

算例

可取字符串

含义

'twister'

'v5uniform'

'v5normal'

'v4'

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

〖说明〗

【例4.3-4】

1）

rngdefault