10349数学建模N6CH04.docx
- 文档编号:23902327
- 上传时间:2023-05-22
- 格式:DOCX
- 页数:31
- 大小:214.46KB
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=
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 10349 数学 建模 N6CH04