matlab课后习题答案第四章.docx

资源描述

matlab课后习题答案第四章.docx

《matlab课后习题答案第四章.docx》由会员分享，可在线阅读，更多相关《matlab课后习题答案第四章.docx（40页珍藏版）》请在冰点文库上搜索。

matlab课后习题答案第四章.docx

matlab课后习题答案第四章

--本页仅作为文档封面，使用时请直接删除即可----内页可以根据需求调整合适字体及大小--

matlab课后习题答案第四章（总22页）

第4章数值运算

习题4及解答

11根据题给的模拟实际测量数据的一组

和

试用数值差分diff或数值梯度gradient指令计算

，然后把

和

曲线绘制在同一张图上，观察数值求导的后果。

（模拟数据从获得）

〖目的〗

强调：

要非常慎用数值导数计算。

练习mat数据文件中数据的获取。

实验数据求导的后果

把两条曲线绘制在同一图上的一种方法。

〖解答〗

（1）从数据文件获得数据的指令

假如文件在当前目录或搜索路径上

clear

load

（2）用diff求导的指令

dt=t

（2）-t

（1）;

yc=diff（y）/dt;%注意yc的长度将比y短1

plot（t,y,'b',t（2:

end）,yc,'r'）

gridon

（3）用gradent求导的指令（图形与上相似）

dt=t

（2）-t

（1）;

yc=gradient（y）/dt;

plot（t,y,'b',t,yc,'r'）

gridon

〖说明〗

不到万不得已，不要进行数值求导。

假若一定要计算数值导数，自变量增量dt要取得比原有数据相对误差高1、2个量级以上。

求导会使数据中原有的噪声放大。

12采用数值计算方法，画出

在

区间曲线，并计算

。

〖提示〗

指定区间内的积分函数可用cumtrapz指令给出。

在计算要求不太高的地方可用find指令算得。

〖目的〗

指定区间内的积分函数的数值计算法和cumtrapz指令。

find指令的应用。

〖解答〗

dt=1e-4;

t=0:

dt:

10;

t=t+（t==0）*eps;

f=sin（t）./t;

s=cumtrapz（f）*dt;

plot（t,s,'LineWidth',3）

ii=find（t==;

s45=s（ii）

s45=

13求函数

的数值积分

，并请采用符号计算尝试复算。

〖提示〗

数值积分均可尝试。

符号积分的局限性。

〖目的〗

符号积分的局限性。

〖解答〗

dx=pi/2000;

x=0:

dx:

pi;

s=trapz（exp（sin（x）.^3））*dx

符号复算的尝试

symsx

f=exp（sin（x）^3）;

ss=int（f,x,0,pi）

Warning:

Explicitintegralcouldnotbefound.

>Inat58

ss=

int（exp（sin（x）^3）,x=0..pi）

14用quad求取

的数值积分，并保证积分的绝对精度为

。

〖目的〗

quadl，精度可控，计算较快。

近似积分指令trapz获得高精度积分的内存和时间代价较高。

〖解答〗

%精度可控的数值积分

fx=@（x）exp（-abs（x））.*abs（sin（x））;

formatlong

sq=quadl（fx,-10*pi,*pi,1e-7）

sq=

%近似积分算法

x=linspace（-10*pi,*pi,1e7）;

dx=x

（2）-x

（1）;

st=trapz（exp（-abs（x））.*abs（sin（x）））*dx

st=

%符号积分算法

y='exp（-abs（x））*abs（sin（x））'

si=vpa（int（y,-10*pi,*pi）,16）

exp（-abs（x））*abs（sin（x））

si=

15求函数

在区间

中的最小值点。

〖目的〗

理解极值概念的邻域性。

如何求最小值。

学习运用作图法求极值或最小值。

感受符号法的局限性。

〖解答〗

（1）采用fminbnd找极小值点

在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。

然后从一系列极小值点中，确定最小值点。

clear

ft=@（t）sin（5*t）.^2.*exp*t.*t）+*abs（t+*t.*cos（2*t）;

disp（'计算中，把[-5,5]分成若干搜索子区间。

'）

N=input（'请输入子区间数N，注意使N>=1'）;%该指令只能在指令窗中运行

tt=linspace（-5,5,N+1）;

fork=1:

[tmin（k）,fobj（k）]=fminbnd（ft,tt（k）,tt（k+1））;

end

[fobj,ii]=sort（fobj）;%将目标值由小到大排列

tmin=tmin（ii）;%使极小值点做与目标值相应的重新排列

fobj,tmin

（2）最后确定的最小值点

在

的不同分割下，经观察，最后确定出

最小值点是

相应目标值是

（3）采用作图法近似确定最小值点（另一方法）

（A）在指令窗中运行以下指令：

clear

ft=@（t）sin（5*t）.^2.*exp*t.*t）+*abs（t+*t.*cos（2*t）;

t=-5:

ff=ft（t）;

plot（t,ff）

gridon,shg

（B）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据

[tmin2,fobj2]=ginput

（1）

tmin2=

fobj2=

出现具有相同数值的刻度区域表明已达最小可分辨状态

（4）符号法求最小值的尝试

symst

fts=sin（5*t）^2*exp*t*t）*t*cos（2*t）+*abs（t+;

dfdt=diff（fts,t）;%求导函数

tmin=solve（dfdt,t）%求导函数的零点

fobj3=subs（fts,t,tmin）%得到一个具体的极值点

tmin=

fobj3=

.024

〖说明〗

最小值是对整个区间而言的，极小值是对邻域而言的。

在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。

这样可以避免把极小值点误作为最小值点。

最小值点是从一系列极小值点和边界点的比较中确定的。

作图法求最小值点，很直观。

假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。

符号法在本例中，只求出一个极值点。

其余很多极值点无法秋初，更不可能得到最小值。

16设

，用数值法和符号法求

。

〖目的〗

学习如何把高阶微分方程写成一阶微分方程组。

ode45解算器的导数函数如何采用匿名函数形式构成。

如何从ode45一组数值解点，求指定自变量对应的函数值。

〖解答〗

（1）改写高阶微分方程为一阶微分方程组

令

，于是据高阶微分方程可写出

（2）运行以下指令求y（t）的数值解

formatlong

ts=[0,1];

y0=[1;0];

dydt=@（t,y）[y

（2）;-2*y

（1）+3*y

（2）+1];%<4>

%匿名函数写成的ode45所需得导数函数

[tt,yy]=ode45（dydt,ts,y0）;

y_05=interp1（tt,yy（:

1）,,'spline'）,%用一维插值求y

y_05=

（3）符号法求解

symst;

ys=dsolve（'D2y-3*Dy+2*y=1','y（0）=1,Dy（0）=0','t'）

ys_05=subs（ys,t,sym（''））

ys=

1/2-1/2*exp（2*t）+exp（t）

ys_05=

.290

〖说明〗

第<4>条指令中的导数函数也可采用M函数文件表达，具体如下。

functionS=prob_DyDt（t,y）

S=[y

（2）;-2*y

（1）+3*y

（2）+1];

17已知矩阵A=magic（8），

（1）求该矩阵的“值空间基阵”B；

（2）写出“A的任何列可用基向量线性表出”的验证程序（提示：

利用rref检验）。

〖目的〗

体验矩阵值空间的基向量组的不唯一性，但它们可以互为线性表出。

利用rref检验两个矩阵能否互为表出。

〖解答〗

（1）A的值空间的三组不同“基”

A=magic（8）;%采用8阶魔方阵作为实验矩阵

[R,ci]=rref（A）;

B1=A（:

ci）%直接从A中取基向量

B2=orth（A）%求A值空间的正交基

[V,D]=eig（A）;

rv=sum（sum（abs（D））>1000*eps）;%非零特征值数就是矩阵的秩

B3=V（:

rv）%取A的非零特征值对应的特征向量作基

B1=

6423

95554

174746

402627

323435

412322

491514

85859

B2=

B3=

（2）验证A的任何列可用B1线性表出

B1_A=rref（[B1,A]）%若B1_A矩阵的下5行全为0，

%就表明A可以被B1的3根基向量线性表出

B1_A=

10010011001

01001034-3-47

001001-3-445-7

00000000000

B2_A=rref（[B2,A]）

B2_A=

Columns1through7

0000000

Columns8through11

0000

B3_A=rref（[B3,A]）

B3_A=

Columns1through7

0000000

Columns8through11

0000

〖说明〗

magic（n）产生魔方阵。

魔方阵具有很多特异的性质。

就其秩而言，当n为奇数时，该矩阵满秩；当n为4的倍数时，矩阵的秩总是3；当n为偶数但不是4倍数时，则矩阵的秩等于（n/2+2）。

关于魔方阵的有关历史，请见第节。

18已知由MATLAB指令创建的矩阵A=gallery（5），试对该矩阵进行特征值分解，并通过验算观察发生的现象。

〖目的〗

展示特征值分解可能存在的数值问题。

condeig是比较严谨的特征值分解指令。

Jordan分解的作用。

〖解答〗

（1）特征值分解

A=gallery（5）

[V,D]=eig（A）;

diag（D）'%为紧凑地显示特征值而写

-911-2163-252

70-69141-4211684

-575575-11493451-13801

3891-38917782-2334593365

1024-10242048-614424572

ans=

Columns1through4

-+-

Column5

（2）验算表明相对误差较大

AE=V*D/V

er_AE=norm（A-AE,'fro'）/norm（A,'fro'）%相对F范数

AE=

+004*

Columns1through4

+-+-

-+-+

+-+-

-+-+

Column5

er_AE=

（3）一个更严谨的特征值分解指令

[Vc,Dc,eigc]=condeig（A）%eigc中的高值时，说明相应的特征值不可信。

Vc=

Columns1through4

+-+

-+-

Column5

Dc=

Columns1through4

000

0+00

00-0

000+

0000

Column5

eigc=

+011*

（4）对A采用Jordan分解并检验

[VJ,DJ]=jordan（A）;%求出准确的广义特征值，使A*VJ=VJ*D成立。

AJ=VJ*DJ/VJ

er_AJ=norm（A-AJ,'fro'）/norm（A,'fro'）

DJ=

01000

00100

00010

00001

00000

AJ=

+004*

er_AJ=

〖说明〗

指令condeig的第3输出量eigc给出的是所谓的“矩阵特征值条件数”。

当特征条件数与

相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大于“几何重数”的特征值。

此时，实施Jordan分解更适宜。

顺便指出：

借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不同概念。

前者描述特征值的问题，后者描述矩阵逆的问题。

本例矩阵A的特征值条件数很高，表明分解不可信。

验算也表明，相对误差较大。

当对矩阵A进行Jordan分解时，可以看到，A具有5重根。

当对Jordan分解进行验算时，相对误差很小。

19求矩阵

的解，A为3阶魔方阵，b是

的全1列向量。

〖提示〗

了解magic指令

rref用于方程求解。

矩阵除法和逆阵法解方程。

〖目的〗

满秩方阵求解的一般过程。

rref用于方程求解。

矩阵除法和逆阵法解方程。

〖解答〗

A=magic（3）;%产生3阶魔方阵

b=ones（3,1）;%（3*1）全1列向量

[R,C]=rref（[A,b]）%GaussJordan消去法解方程，同时判断解的唯一性

x=A\b%矩阵除解方程

xx=inv（A）*b%逆阵法解方程

123

xx=

〖说明〗

rref指令通过对增广矩阵进行消去法操作完成解方程。

由分解得到的3根“坐标向量”和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。

在本例情况下，矩阵除、逆阵法、rref法所得解一致。

110求矩阵

的解，A为4阶魔方阵，b是

的全1列向量。

〖提示〗

用rref可观察A不满秩，但b在A的值空间中，这类方程用无数解。

矩阵除法能正确求得这类方程的特解。

逆阵法不能求得这类方程的特解。

注意特解和齐次解

〖目的〗

A不满秩，但b在A的值空间中，这类方程的求解过程。

rref用于方程求解。

矩阵除法能正确求得这类方程的特解。

逆阵法不能求得这类方程的特解。

解的验证方法。

齐次解的获取。

全解的获得。

〖解答〗

（1）借助增广矩阵用指令rref求解

A=magic（4）;%产生3阶魔方阵

b=ones（4,1）;%全1列向量

[R,C]=rref（[A,b]）%求解，并判断解的唯一性

00000

123

关于以上结果的说明：

R阶梯阵提供的信息

前4列是原A阵经消元变换后的阶梯阵；而第5列是原b向量经相同变换后的结果。

R的前三列为“基”，说明原A阵秩为3；而第4列的前三个元素，表示原A阵的第4列由其前三列线性组合而成时的加权系数，即方程的一个解。

R的第5列表明：

b可由原A阵的前三列线性表出；b给出了方程的一个解；由于原A阵“缺秩”，所以方程的确解不唯一。

C数组提供的信息

该数组中的三个元素表示变换取原A阵的第1，2，3列为基。

该数组的元素总数就是“原A阵的秩”

（2）矩阵除求得的解

x=A\b

Warning:

Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=.

运行结果指示：

矩阵除法给出的解与rref解相同。

（实际上，MATLAB在设计“除法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。

）

（3）逆阵法的解

xx=inv（A）*b

Warning:

Matrixisclosetosingularorbadlyscaled.

Resultsmaybeinaccurate.RCOND=.

xx=

（3）验证前面所得的解

b_rref=A（:

C）*R（C,5）%验算rref的解

b_d=A*x%验算矩阵除的解

b_inv=A*xx%验算逆阵法的解

b_rref=

b_d=

b_inv=

显然，在本例中，逆阵法的解是错误的。

原因是：

A不满秩，A的逆阵在理论上不存在。

这里所给出的逆阵是不可信的。

（4）求齐次解

xg=null（A）%Ax=0的齐次解

xg=

（5）方程的全解

齐次解的任何倍与特解之和就构成方程的全解。

下面通过一组随机系数验证。

rngdefault%为本书结果可被读者核对而设，并非必要。

f=randn（1,6）%6个随机系数

xx=repmat（x,1,6）+xg*f%产生6个不同的特解

A*xx%所得结果的每列都应该是全1，即等于b.

xx=

ans=

〖解答〗

（在用除法和逆阵法求解时出现）警告信息中RCOND=是矩阵A的估计条件倒数。

该数愈接近0，A就愈“病态”；该数接近1时，A就愈“良态”。

该条件数由rcond（A）给出。

注意：

rcond条件倒数与cond条件数的算法不同。

111求矩阵

的解，A为4阶魔方阵，

。

〖提示〗

由rref可以看出A不满秩，b不在A的值空间中，方程没有准确解。

但可求最小二乘近似解。

〖目的〗

A不满秩，b不在A的值空间中，方程没有准确解。

〖解答〗

（1）借助增广矩阵用指令rref求解

A=magic（4）;%产生3阶魔方阵

b=（1:

4）';

[R,C]=rref（[A,b]）%求解，并判断解的唯一性

10010

01030

001-30

00001

1235

（2）用伪逆求最小二乘近似解（超出范围，仅供参考。

）

x=pinv（A）*b%非准确解

b_pinv=A*x%验算

b_pinv=

〖解答〗

C表明，A的秩为3，A不满秩；R第5列第4元素非零，说明b不在A的值空间中，因此方程没有准确解，但可以求最小二乘近似解。

112求

的实数解。

〖提示〗

在适当范围内，作图观察一元复杂函数的形态：

观察解的存在性；解的唯一性。

进而，借助图形法求近似解。

匿名函数的使用方法。

fzero指令的用法。

〖目的〗

作图法求一元复杂函数解上的作用：

观察解的存在性；解的唯一性；得近似解。

匿名函数的使用方法。

fzero指令的用法。

〖解答〗

（1）作图观察函数并求近似解

t=-1:

y=@（t）+t-10*exp*t）.*abs（sin（sin（t）））;

plot（t,y（t））%利用匿名函数求y函数值

gridon,shg

[tt1,yy1]=ginput

（1）%从图形获得近似解

tt1=

yy1=

（2）进一步利用fzero求精确解

[t,yt]=fzero（y,tt1）

yt=

〖说明〗

假如在从图上获取数据前，先把零点附近图形放大，可以得到精度更高的近似解。

113求解二元函数方程组

的解。

〖目的〗

solve指令的用法。

〖解答〗

（1）符号法只能得到两组解

S=solve（'sin（x-y）','cos（x+y）','x','y'）;

X=,Y=

[1/4*pi]

[-1/4*pi]

[1/4*pi]

[-1/4*pi]

（2）数值法可以看到许多解，并从中找到规律

x=-6:

6;y=x;[X,Y]=meshgrid（x,y）;

Z1=sin（X-Y）;

Z2=cos（X+Y）;

contour（X,Y,Z1,3）%在Z1的取值范围内取3个等分点作为等位线取值。

%由于Z1关于z1=0对称，所以中间线，即Z1=0的等位线。

axissquare

holdon

contour（X,Y,Z2,3）%注意：

所有绿线交点都是解

holdoff

gridon;shg

114假定某窑工艺瓷器的烧制成品合格率为，现该窑烧制100件瓷器，请画出合格产品数的概率分布曲线。

〖目的〗

指令binopdf的用法。

绘图指令stem的用法。

〖解答〗

N=100;p=;%给定二项分布的特征参数

k=0:

N;%定义事件A发生的次数数组

pdf=binopdf（k,N,p）;%算出各发生次数的概率

stem（k,pdf）

gridon

shg

展开阅读全文