10349数学建模N6CH04.docx
《10349数学建模N6CH04.docx》由会员分享，可在线阅读，更多相关《10349数学建模N6CH04.docx（31页珍藏版）》请在冰豆网上搜索。
10349数学建模N6CH04
第四章数值计算
.1线性方程组的解
.1.1LU分解、行列式、逆和恰定方程的解
【例4.1-1】“求逆”法和“左除”法解恰定方程的性能对比
（1）
randn（'state',0）;
A=gallery（'randsvd',100,2e13,2）;
x=ones（100,1）;
b=A*x;
cond（A）
ans=
1.9990e+013
（2）
tic
xi=inv（A）*b;
ti=toc
eri=norm（x-xi）
rei=norm（A*xi-b）/norm（b）
ti=
0.4400
eri=
0.0469
rei=
0.0047
（3）
tic;xd=A\b;
td=toc,
erd=norm（x-xd）,
red=norm（A*xd-b）/norm（b）
td=
0.0600
erd=
0.0078
red=
2.6829e-015
.1.2奇异值分解和矩阵结构
.1.3线性二乘问题的解
【例4.1-2】对于超定方程
，进行三种解法比较。
其中
取MATLAB库中的特殊函数生成。
A=gallery（5）;A（:
1）=[];y=[1.77.56.30.83-0.082]';
x=inv（A'*A）*A'*y,xx=pinv（A）*y,xxx=A\y
Warning:
Matrixisclosetosingularorbadlyscaled.
Resultsmaybeinaccurate.RCOND=7.087751e-018.
x=
3.4811
5.1595
0.9534
-0.0466
xx=
3.4759
5.1948
0.7121
-0.1101
Warning:
Rankdeficient,rank=3tol=1.0829e-010.
xxx=
3.4605
5.2987
0
-0.2974
nx=norm（x）,nxx=norm（xx）,nxxx=norm（xxx）
nx=
6.2968e+000
nxx=
6.2918e+000
nxxx=
6.3356e+000
e=norm（y-A*x）,ee=norm（y-A*xx）,eee=norm（y-A*xxx）
e=
1.1020e+000
ee=
4.7424e-002
eee=
4.7424e-002
.2特征值分解和矩阵函数
.2.1特征值分解问题
【例4.2-1】简单实阵的特征值问题。
A=[1,-3;2,2/3];[V,D]=eig（A）
V=
-0.7728+0.0527i-0.7728-0.0527i
0+0.6325i0-0.6325i
D=
0.8333+2.4438i0
00.8333-2.4438i
【例4.2-2】把例4.2-1中的复数特征值对角阵D转换成实数块对角阵，使VR*DR/VR=A。
[VR,DR]=cdf2rdf（V,D）
VR=
-0.77280.0527
00.6325
DR=
0.83332.4438
-2.44380.8333
【例4.2-3】对亏损矩阵进行Jordan分解。
A=gallery（5）
[VJ,DJ]=jordan（A）;
[V,D,c_eig]=condeig（A）;c_equ=cond（A）;
DJ,D,c_eig,c_equ
A=
-911-2163-252
70-69141-4211684
-575575-11493451-13801
3891-38917782-2334593365
1024-10242048-614424572
DJ=
01000
00100
00010
00001
00000
D=
Columns1through4
-0.0328+0.0243i000
0-0.0328-0.0243i00
000.0130+0.0379i0
0000.0130-0.0379i
0000
Column5
0
0
0
0
0.0396
c_eig=
1.0e+010*
2.1016
2.1016
2.0251
2.0251
1.9796
c_equ=
5.2129e+017
.2.2矩阵的谱分解和矩阵函数
【例4.2-4】数组乘方与矩阵乘方的比较。
clear,A=[123;456;789];
A_Ap=A.^0.3
A_Mp=A^0.3
A_Ap=
1.00001.23111.3904
1.51571.62071.7118
1.79281.86611.9332
A_Mp=
0.6962+0.6032i0.4358+0.1636i0.1755-0.2759i
0.6325+0.0666i0.7309+0.0181i0.8292-0.0305i
0.5688-0.4700i1.0259-0.1275i1.4830+0.2150i
【例4.2-5】标量的数组乘方和矩阵乘方的比较。
（A取自例4.2-4）
pA_A=（0.3）.^A
pA_M=（0.3）^A
pA_A=
0.30000.09000.0270
0.00810.00240.0007
0.00020.00010.0000
pA_M=
2.93420.4175-1.0993
-0.02780.7495-0.4731
-1.9898-0.91841.1531
【例4.2-6】sin的数组运算和矩阵运算比较。
（A取自例4.2-4）
A_sinA=sin（A）
A_sinM=funm（A,'sin'）
A_sinA=
0.84150.90930.1411
-0.7568-0.9589-0.2794
0.65700.98940.4121
A_sinM=
-0.6928-0.23060.2316
-0.1724-0.1434-0.1143
0.3479-0.0561-0.4602
.3多项式和卷积
.3.1多项式
10一多项式表达方式的约定
10二多项式运算函数
【例4.3-1】求
的“商”及“余”多项式。
p1=conv（[1,0,2],conv（[1,4],[1,1]））;
p2=[1011];
[q,r]=deconv（p1,p2）;
cq='商多项式为';cr='余多项式为';
disp（[cq,poly2str（q,'s'）]）,disp（[cr,poly2str（r,'s'）]）
商多项式为s+5
余多项式为5s^2+4s+3
【例4.3-2】求3阶方阵A的特征多项式。
A=[111213;141516;171819];
PA=poly（A）
PPA=poly2str（PA,'s'）
PA=
1.0000-45.0000-18.0000-0.0000
PPA=
s^3-45s^2-18s-2.8387e-015
【例4.3.1.2-3】由给定根向量求多项式系数向量。
R=[-0.5,-0.3+0.4*i,-0.3-0.4*i];
P=poly（R）
PR=real（P）
PPR=poly2str（PR,'x'）
P=
1.00001.10000.55000.1250
PR=
1.00001.10000.55000.1250
PPR=
x^3+1.1x^2+0.55x+0.125
【例4.3-4】两种多项式求值指令的差别。
S=pascal（4）
P=poly（S）;PP=poly2str（P,'s'）
PA=polyval（P,S）
PM=polyvalm（P,S）
S=
1111
1234
13610
141020
PP=
s^4-29s^3+72s^2-29s+1
PA=
1.0e+004*
0.00160.00160.00160.0016
0.00160.0015-0.0140-0.0563
0.0016-0.0140-0.2549-1.2089
0.0016-0.0563-1.2089-4.3779
PM=
1.0e-011*
-0.00770.0053-0.00960.0430
-0.00680.0481-0.01100.1222
0.00750.1400-0.00950.2608
0.04300.2920-0.00070.4737
【例4.3-5】部分分式展开。
a=[1,3,4,2,7,2];
b=[3,2,5,4,6];
[r,s,k]=residue（b,a）
r=
1.1274+1.1513i
1.1274-1.1513i
-0.0232-0.0722i
-0.0232+0.0722i
0.7916
s=
-1.7680+1.2673i
-1.7680-1.2673i
0.4176+1.1130i
0.4176-1.1130i
-0.2991
k=
[]
10三拟合和插值
【例4.3-6】对于给定数据对x0,y0，求拟合三阶多项式，并图示拟合情况。
（见图4.3-1）
x0=0:
0.1:
1;
y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
n=3;
P=polyfit（x0,y0,n）
xx=0:
0.01:
1;
yy=polyval（P,xx）;
plot（xx,yy,'-b',x0,y0,'.r','MarkerSize',20）,xlabel（'x'）
P=
56.6915-87.117440.0070-0.9043
图4.3-1采用三次多项式所得的拟合曲线
【例4.3-7】以上例所给数据，研究一维插值，并观察插值与拟合的区别。
x0=0:
0.1:
1;
y0=[-.447,1.978,3.11,5.25,5.02,4.66,4.01,4.58,3.45,5.35,9.22];
xi=0:
0.02:
1;
yi=interp1（x0,y0,xi,'cubic'）;
plot（xi,yi,'-b',x0,y0,'.r','MarkerSize',20）,xlabel（'x'）
图4.3-2通过三次多项式插值所得的曲线
.3.2卷积
10一离散序列的数值卷积
10二MATLAB的“卷积”指令
【例4.3-8】有序列
和
。
（A）求这两个完整序列的卷积，并图示。
（B）假设
序列中最后4个非零值未知，而成为截尾序列，求卷积并图示。
（见图4.3-3）
a=ones（1,10）;n1=3;n2=12;
b=ones（1,8）;n3=2;n4=9;
c=conv（a,b）;nc1=n1+n3;nc2=n2+n4;
kc=nc1:
nc2;
aa=a（1:
6）;nn1=3;nn2=8;
cc=conv（aa,b）;ncc1=nn1+n3;
nx=nn2+n4;
ncc2=min（nn1+n4,nn2+n3）;
kx=（ncc2+1）:
nx;kcc=ncc1:
ncc2;N=length（kcc）;
stem（kcc,cc（1:
N）,'r','filled'）
axis（[nc1-2,nc2+2,0,10]）,grid,holdon
stem（kc,c,'b'）,stem（kx,cc（N+1:
end）,'g'）,holdoff
图4.3-3“完整”序列卷积和“截尾”序列卷积
.4数据分析函数
.4.1随机数发生器和统计分析指令
【例4.4-1】基本统计示例。
randn（'state',0）
A=randn（1000,4）;
AMAX=max（A）,AMIN=min（A）
AMED=median（A）
AMEAN=mean（A）
ASTD=std（A）
AMAX=
2.73163.20253.41283.0868
AMIN=
-2.6442-3.0737-3.5027-3.0461
AMED=
-0.01310.05960.01220.0459
AMEAN=
-0.04310.04550.01770.0263
ASTD=
0.94351.03131.02480.9913
【例4.4-2】cov和corrcoef的使用示例。
rand（'state',1）
X=rand（10,3）;Y=rand（10,3）;
mx=mean（X）;Xmx=X-ones（size（X））*diag（mx）;
CCX=Xmx'*Xmx/（size（Xmx,1）-1）
CX=cov（X）,CY=cov（Y）
Cxy=cov（X,Y）
PX=corrcoef（X）
Pxy=corrcoef（X,Y）
CCX=
0.08710.0242-0.0073
0.02420.08460.0056
-0.00730.00560.0607
CX=
0.08710.0242-0.0073
0.02420.08460.0056
-0.00730.00560.0607
CY=
0.07210.00130.0165
0.00130.0714-0.0535
0.0165-0.05350.0720
Cxy=
0.0761-0.0012
-0.00120.0708
PX=
1.00000.2819-0.1005
0.28191.00000.0785
-0.10050.07851.0000
Pxy=
1.0000-0.0168
-0.01681.0000
.4.2差分和累计指令
【例4.4-3】用一个简单矩阵表现diff和gradient指令计算方式。
F=[1,2,3;4,5,6;7,8,9]
Dx=diff（F）
Dx_2=diff（F,1,2）
[FX,FY]=gradient（F）
[FX_2,FY_2]=gradient（F,0.5）
F=
123
456
789
Dx=
333
333
Dx_2=
11
11
11
FX=
111
111
111
FY=
333
333
333
FX_2=
222
222
222
FY_2=
666
666
666
【例4.4-4】函数
的梯度
和
，用数值计算验证，并图示。
（见图4.4-1）
clear,
dd=0.2;
x=-1:
dd:
1;y=x;
[X,Y]=meshgrid（x,y）;
Z=（X.^2）+（Y.^2）;
[DZx,DZy]=gradient（Z,dd,dd）;
DZ2=4*del2（Z,dd,dd）;
DDZx0=DZx-2*X;DDZy0=DZy-2*Y;DDZ20=DZ2-4;
subplot（1,3,1）,stem3（X,Y,DDZx0）
subplot（1,3,2）,stem3（X,Y,DDZy0）
subplot（1,3,3）,stem3（X,Y,DDZ20）
axis（[-1,1,-1,1,-0.4,0.4]）
xlabel（'x'）,ylabel（'y'）
图4.4-1理论计算和数值计算的差别图示
【例4.4-5】求积分
，其中
，
。
dt=0.1;t=（0:
dt:
10）';
y=exp（-0.8*t.*abs（sin（t）））;
ss=dt*cumsum（y）;
ss10=dt*sum（y）;ssend=ss（end）;
st=cumtrapz（t,y）;
st10=trapz（t,y）;stend=st（end）;
disp（[blanks（5）,'sum',blanks（6）,'cumsum',blanks（4）,'trapz',blanks（5）,'cumtrapz']）
disp（[ss10,ssend,st10,stend]）
plot（t,y,'b:
',t,ss,'r-',t,st,'r.'）
legend（'y（x）','cunsum','cumtrapz',0）
sumcumsumtrapzcumtrapz
2.70822.70822.65762.6576
图4.4-2矩形法和梯形法求积比较
.5MATLAB泛函指令
.5.1求函数零点
【例4.5-1】通过求
的零点，综合叙述相关指令的用法。
P1=0.1;P2=0.5;
y_C='sin（x）.^2.*exp（-P1*x）-P2*abs（x）';
x=-10:
0.01:
10;
Y=eval（y_C）;
clf,plot（x,Y,'r'）;holdon,plot（x,zeros（size（x））,'k'）;
xlabel（'t'）;ylabel（'y（t）'）,holdoff
图4.5-1函数零点分布观察图
zoomon
[tt,yy]=ginput（5）;zoomoff
图4.5-2局部放大和利用鼠标取值图
tt
tt=
-2.0046
-0.5300
-0.0115
0.6106
1.6590
[t4,y4,exitflag]=fzero（y_C,tt（4）,[],P1,P2）
Zerofoundintheinterval:
[0.59333,0.6106].
t4=
0.5993
y4=
0
exitflag=
1
.5.2求函数极值点
【例4.5-2】求
的极小值点。
它即是著名的Rosenbrock's"Banana"测试函数，它的理论极小值是
。
该测试函数有一片浅谷，许多算法难以越过此谷。
ff=inline（'100*（x
（2）-x
（1）^2）^2+（1-x
（1））^2','x'）;
x0=[-1.2,1];[sx,sfval,sexit,soutput]=fminsearch（ff,x0）
Optimizationterminatedsuccessfully:
thecurrentxsatisfiestheterminationcriteriausingOPTIONS.TolXof1.000000e-004
andF（X）satisfiestheconvergencecriteriausingOPTIONS.TolFunof1.000000e-004
sx=
1.00001.0000
sfval=
8.1777e-010
sexit=
1
soutput=
iterations:
85
funcCount:
159
algorithm:
'Nelder-Meadsimplexdirectsearch'
.5.3数值积分
【例4.5-3】求
，其精确值为
。
（1）
[fun.m]
functiony=fun（x）
y=exp（-x.*x）;
（2）
Hf=@fun;
Isim=quad（Hf,0,1）,
IL=quadl（Hf,0,1）
Isim=
0.7468
IL=
0.7468
.5.4解常微分方程
【例4.5-4】求微分方程
在初始条件
情况下的解，并图示。
（见图4.5-3和4.5-4）
（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）;%<3>
plot（tt,yy（:
1））,title（'x（t）'）
图4.5-3微分方程解
（4）
plot（yy（:
1）,yy（:
2））
图4.5-4相平面轨迹
.6信号处理
.6.1快速Fourier变换和逆变换
【例4.6-1】利用fft和ifft指令重新计算两序列的卷积。
所给已知序列为
，
。
clear
a=ones（1,13）;a（[1,2,3]）=0;
b=ones（1,10）;b（[1,2]）=0;
c=conv（a,b）;
M=32;
AF=fft（a,M）;BF=fft（b,M）;
CF=AF.*BF;
cc=real（ifft（CF））;
nn=0:
（M-1）;
c（M）=0;
error=c-cc;
subplot（2,1,1）,stem（nn,c,'fill'）,grid,axis（[0,31,0,9]）
xlabel（'nn'）,ylabel（'cc'）
subplot（2,1,2）,stem（nn,error,'fill'）,axis（[0,31,-1,1]）
ylabel（'error'）
图4.6-1变换法和直接法求卷积结果比较
【例4.6-2】fft在信号分析中的应用。
试用频谱分析方法从受噪声污染的信号
中鉴别出有用信号。
在此，
。
clear,randn（'state',0）
t=linspace（0,10,512）;
y=3*sin（5*t）-6*cos（9*t）+5*randn（size（t））;
plot（t,y）
图4.6-2受噪声污染的信号
Y=
升级会员 