当前位置：搜档网 › matlab课后习题答案第四章

matlab课后习题答案第四章

第4章数值运算

习题 4 及解答

1 根据题给的模拟实际测量数据的一组t和)(t y试用数值差分diff

或数值梯度gradient指令计算)(t

y'曲线绘制

y'，然后把)(t y和)(t 在同一张图上，观察数值求导的后果。（模拟数据从prob_data401.mat获得）

〖目的〗

●强调：要非常慎用数值导数计算。

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

●实验数据求导的后果

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

〖解答〗

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

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

clear

load prob_data401.mat

（2）用diff求导的指令

dt=t(2)-t(1);

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

plot(t,y,'b',t(2:end),yc,'r')

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

dt=t(2)-t(1);

yc=gradient(y)/dt;

plot(t,y,'b',t,yc,'r') grid on

〖说明〗

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

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

以上。

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

2 采用数值计算方法，画出dt t

x y x

sin )(在]10 ,0[区间曲线，并计算)5.4(y 。

〖提示〗

● 指定区间内的积分函数可用cumtrapz 指令给出。 ● )5.4(y 在计算要求不太高的地方可用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==4.5); s45=s(ii) s45 = 1.6541

3 求函数x

x f 3sin )(=的数值积分?=π

0 )(dx x f s ，并请采用符号计算

尝试复算。

〖提示〗

● 数值积分均可尝试。 ● 符号积分的局限性。

〖目的〗

● 符号积分的局限性。〖解答〗

dx=pi/2000; x=0:dx:pi;

s=trapz(exp(sin(x).^3))*dx s =

5.1370

符号复算的尝试

syms x

f=exp(sin(x)^3); ss=int(f,x,0,pi)

Warning: Explicit integral could not be found. > In sym.int at 58 ss =

int(exp(sin(x)^3),x = 0 .. pi)

4 用

quad 求取dx x e x sin 7.15?--π

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

度

为

10-。

〖目的〗

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

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

%精度可控的数值积分

fx=@(x)exp(-abs(x)).*abs(sin(x)); format long

sq=quadl(fx,-10*pi,1.7*pi,1e-7) sq =

1.08784993815498

%近似积分算法

x=linspace(-10*pi,1.7*pi,1e7); dx=x(2)-x(1);

st=trapz(exp(-abs(x)).*abs(sin(x)))*dx st =

1.08784949973430

%符号积分算法

y='exp(-abs(x))*abs(sin(x))'

si=vpa(int(y,-10*pi,1.7*pi),16) y =

exp(-abs(x))*abs(sin(x)) si =

1.087849499412911

5 求函数5.08.12cos 5.1)

5(sin )(2

06.02

++-=t t t e t t f t 在区间]5,5[-中的最小

值点。

〖目的〗

● 理解极值概念的邻域性。 ● 如何求最小值。

● 学习运用作图法求极值或最小值。 ● 感受符号法的局限性。〖解答〗

（1）采用fminbnd 找极小值点在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。然后从一系列极小值点中，确定最小值点。 clear

ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t); disp('计算中，把[ -5,5] 分成若干搜索子区间。')

N=input(' 请输入子区间数 N ，注意使N>=1 ？');%该指令只能在指令窗中运行 t t=linspace(-5,5,N+1); f or k=1:N

[tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1));

e nd [fobj,ii]=sort(fobj); %将目标值由小到大排列 t min=tmin(ii); %使极小值点做与目标值相应的重新排列 fobj,tmin

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

在10,,2,1 N 的不同分割下，经观察，最后确定出最小值点是 -1.28498111480531 相应目标值是 -0.18604801006545

（3）采用作图法近似确定最小值点（另一方法）（A ）在指令窗中运行以下指令： clear

ft=@(t)sin(5*t).^2.*exp(0.06*t.*t)+1.8*abs(t+0.5)-1.5*t.*cos(2*t); t=-5:0.001:5; ff=ft(t); plot(t,ff) grid on,shg （B ）经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线”，利用它点击极值点可得到最小值数据

[tmin2,fobj2]=ginput(1) tmin2 =

-1.28500000993975 fobj2 =

-0.18604799369136

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

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

syms t

fts=sin(5*t)^2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8*abs(t+0.5);

dfdt=diff(fts,t); %求导函数

tmin=solve(dfdt,t) %求导函数的零点

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

tmin =

-.60100931947716486053884417850955e-2

fobj3 =

.89909908144684551670208797723124

〖说明〗

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

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

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

●作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小

值点有相当好的精度。

●符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。

6 设0)0(,1)0(,1)(2)(3)

(22===+-dt dy y t y dt t dy dt t y d ，用数值法和符号法求

5.0)(=t t y 。

〖目的〗

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

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

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

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

令dt

t dy t y t y t y )

()(),()(21=

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

++-==1)(3)(2)

()()(21221t y t y dt

t dy t y dt t dy

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

format long 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),0.5,'spline'), %用一维插值求y(0.5) y_05 =

0.78958020790127

（3）符号法求解

syms t;

ys=dsolve('D2y-3*Dy+2*y=1','y(0)=1,Dy(0)=0','t') ys_05=subs(ys,t,sym('0.5')) ys =

1/2-1/2*exp(2*t)+exp(t) ys_05 =

.78958035647060552916850705213780

〖说明〗

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

function S=prob_DyDt(t,y) S=[y(2);-2*y(1)+3*y(2)+1];

7 已知矩阵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(:,1:rv) %取A的非零特征值对应的特征向量作基

B1 =

64 2 3

9 55 54

17 47 46

40 26 27

32 34 35

41 23 22

49 15 14

8 58 59

B2 =

-0.3536 0.5401 0.3536

-0.3536 -0.3858 -0.3536

-0.3536 -0.2315 -0.3536

-0.3536 0.0772 0.3536

-0.3536 -0.0772 0.3536

-0.3536 0.2315 -0.3536

-0.3536 0.3858 -0.3536

-0.3536 -0.5401 0.3536

B3 =

0.3536 0.6270 0.3913

0.3536 -0.4815 -0.2458

0.3536 -0.3361 -0.1004

0.3536 0.1906 -0.0451

0.3536 0.0451 -0.1906

0.3536 0.1004 0.3361

0.3536 0.2458 0.4815

0.3536 -0.3913 -0.6270

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

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

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

B1_A =

1 0 0 1 0 0 1 1 0 0 1

0 1 0 0 1 0 3 4 -3 -4 7

0 0 1 0 0 1 -3 -4 4 5 -7

0 0 0 0 0 0 0 0 0 0 0

B2_A=rref([B2,A])

B2_A =

Columns 1 through 7

1.0000 0 0 -91.9239 -91.9239 -91.9239 -91.9239

0 1.0000 0 51.8459 -51.8459 -51.8459 51.8459

0 0 1.0000 9.8995 -7.0711 -4.2426 1.4142

0 0 0 0 0 0 0

Columns 8 through 11

-91.9239 -91.9239 -91.9239 -91.9239

51.8459 -51.8459 -51.8459 51.8459

-1.4142 4.2426 7.0711 -9.8995

0 0 0 0

B3_A=rref([B3,A])

B3_A =

Columns 1 through 7

1.0000 0 0 91.9239 91.9239 91.9239 91.9239

0 1.0000 0 42.3447 -38.1021 -33.8594 29.6168

0 0 1.0000 12.6462 -16.8889 -21.1315 25.3741

0 0 0 0 0 0 0

Columns 8 through 11

91.9239 91.9239 91.9239 91.9239

25.3741 -21.1315 -16.8889 12.6462

29.6168 -33.8594 -38.1021 42.3447

0 0 0 0

〖说明〗

●magic(n)产生魔方阵。魔方阵具有很多特异的性质。就其秩而言，当n为奇数时，该矩

阵满秩；当n 为4的倍数时，矩阵的秩总是3；当 n 为偶数但不是4倍数时，则矩阵的秩等于 (n/2+2)。关于魔方阵的有关历史，请见第6.1.3节。

8 已知由MATLAB指令创建的矩阵A=gallery(5)，试对该矩阵进行

特征值分解，并通过验算观察发生的现象。

〖目的〗

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

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

●Jordan分解的作用。

〖解答〗

（1）特征值分解

A=gallery(5)

[V,D]=eig(A);

diag(D)' %为紧凑地显示特征值而写

A =

-9 11 -21 63 -252

70 -69 141 -421 1684

-575 575 -1149 3451 -13801

3891 -3891 7782 -23345 93365

1024 -1024 2048 -6144 24572

ans =

Columns 1 through 4

-0.0181 -0.0054 - 0.0171i -0.0054 + 0.0171i 0.0144 - 0.0104i

Column 5

0.0144 + 0.0104i

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

AE=V*D/V

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

AE =

1.0e+004 *

Columns 1 through 4

-0.0009 + 0.0000i 0.0011 - 0.0000i -0.0021 + 0.0000i 0.0063 - 0.0000i

0.0070 - 0.0000i -0.0069 + 0.0000i 0.0141 - 0.0000i -0.0421 + 0.0000i

-0.0575 + 0.0000i 0.0575 - 0.0000i -0.1149 + 0.0000i 0.3451 - 0.0000i

0.3891 - 0.0000i -0.3891 + 0.0000i 0.7781 - 0.0000i -2.3343 + 0.0000i

0.1024 - 0.0000i -0.1024 + 0.0000i 0.2048 - 0.0000i -0.6144 + 0.0000i

Column 5

-0.0252 + 0.0000i

0.1684 - 0.0000i

-1.3800 + 0.0000i

9.3359 - 0.0001i

2.4570 - 0.0000i

er_AE =

6.9310e-005

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

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

Vc =

Columns 1 through 4

-0.0000 -0.0000 + 0.0000i -0.0000 - 0.0000i 0.0000 + 0.0000i

0.0206 0.0207 + 0.0000i 0.0207 - 0.0000i 0.0207 + 0.0000i

-0.1397 -0.1397 + 0.0000i -0.1397 - 0.0000i -0.1397 + 0.0000i

0.9574 0.9574 0.9574 0.9574

0.2519 0.2519 - 0.0000i 0.2519 + 0.0000i 0.2519 - 0.0000i

Column 5

0.0000 - 0.0000i

0.0207 - 0.0000i

-0.1397 - 0.0000i

0.9574

0.2519 + 0.0000i

Dc =

Columns 1 through 4

-0.0181 0 0 0

0 -0.0054 + 0.0171i 0 0

0 0 -0.0054 - 0.0171i 0

0 0 0 0.0144 + 0.0104i

0 0 0 0

Column 5

0.0144 - 0.0104i

eigc =

1.0e+011 *

5.2687

5.2313

5.1725

5.1724

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

0 1 0 0 0

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

0 0 0 0 0

AJ =

1.0e+004 *

-0.0009 0.0011 -0.0021 0.0063 -0.0252

0.0070 -0.0069 0.0141 -0.0421 0.1684

-0.0575 0.0575 -0.1149 0.3451 -1.3801

0.3891 -0.3891 0.7782 -2.3345 9.3365

0.1024 -0.1024 0.2048 -0.6144 2.4572

er_AJ =

2.0500e-011

〖说明〗

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

1相当时，就意味着矩阵A可能“退化”，即矩阵可能存在“代数重数”大数与

eps

于“几何重数”的特征值。此时，实施Jordan分解更适宜。

●顺便指出：借助condeig算得的特征值条件数与cond指令算得的矩阵条件数是两个不

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

●本例矩阵A的特征值条件数很高，表明分解不可信。验算也表明，相对误差较大。

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

时，相对误差很小。

9 求矩阵b

3(?的全1列向量。

Ax=的解，A为3阶魔方阵，b是)1

〖提示〗

●了解magic指令

●rref 用于方程求解。

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

〖目的〗

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

●rref 用于方程求解。

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

〖解答〗

A=magic(3); %产生3阶魔方阵

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

[R,C]=rref([A,b]) %Gauss Jordan消去法解方程，同时判断解的唯一性

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

xx=inv(A)*b %逆阵法解方程

R =

1.0000 0 0 0.0667

0 1.0000 0 0.0667

0 0 1.0000 0.0667

C =

1 2 3

x =

0.0667

xx =

0.0667

〖说明〗

●rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的3根“坐标向量”

和（或）C3指示的3根基向量，可见A3满秩，因此方程解唯一。

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

10 求矩阵b

4(?的全1列向量。

Ax=的解，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]) %求解，并判断解的唯一性

R =

1.0000 0 0 1.0000 0.0588

0 1.0000 0 3.0000 0.1176

0 0 1.0000 -3.0000 -0.0588

0 0 0 0 0

C =

1 2 3

关于以上结果的说明：

●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: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 1.306145e-017.

x =

0.0588

0.1176

-0.0588

运行结果指示：矩阵除法给出的解与rref解相同。（实际上，MATLAB在设计“除法”程序时，针对“b在A值空间中”的情况，就是用rref求解的。）

（3）逆阵法的解

xx=inv(A)*b

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 1.306145e-017.

xx =

0.0469

0.1875

-0.0625

-0.0156

（3）验证前面所得的解

b_rref=A(:,C)*R(C,5) %验算rref的解

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

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

b_rref =

b_d =

b_inv =

0.7344

1.1719 1.8594

显然，在本例中，逆阵法的解是错误的。原因是：A 不满秩，A 的逆阵在理论上不存在。这里所给出的逆阵是不可信的。

（4）求齐次解

xg=null(A)

%Ax=0的齐次解

xg =

0.2236 0.6708 -0.6708 -0.2236

（5）方程的全解齐次解的任何倍与特解之和就构成方程的全解。下面通过一组随机系数验证。

rng default %为本书结果可被读者核对而设，并非必要。 f=randn(1,6) %6个随机系数 xx=repmat(x,1,6)+xg*f %产生6个不同的特解 A*xx %所得结果的每列都应该是全1，即等于b. f =

0.5377 1.8339 -2.2588 0.8622 0.3188 -1.3077 xx =

0.1790 0.4689 -0.4463 0.2516 0.1301 -0.2336 0.4783 1.3479 -1.3976 0.6960 0.3315 -0.7596 -0.4195 -1.2890 1.4565 -0.6372 -0.2727 0.8184 -0.1202 -0.4101 0.5051 -0.1928 -0.0713 0.2924 ans =

1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000

〖解答〗

● （在用除法和逆阵法求解时出现）警告信息中RCOND = 1.306145e-017是矩阵A 的估计

条件倒数。该数愈接近0，A 就愈“病态”；该数接近1时，A 就愈“良态”。该条件数由rcond(A)给出。注意：rcond 条件倒数与cond 条件数的算法不同。

11 求矩阵b Ax =的解，A 为4阶魔方阵，?????

???????=4321b 。〖提示〗

● 由rref 可以看出A 不满秩，b 不在A 的值空间中，方程没有准确解。 ● 但可求最小二乘近似

解。

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

〖解答〗

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

A=magic(4); %产生3阶魔方阵

b=(1:4)';

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

R =

1 0 0 1 0

0 1 0 3 0

0 0 1 -3 0

0 0 0 0 1

C =

1 2 3 5

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

x=pinv(A)*b %非准确解

b_pinv=A*x %验算

x =

0.0235

0.1235

0.0235

b_pinv =

1.3000

2.9000

2.1000

3.7000

〖解答〗

●C表明，A的秩为3，A不满秩；R第5列第4元素非零，说明b不在A的值空间中，因

此方程没有准确解，但可以求最小二乘近似解。

12 求0

--t

t t的实数解。

5.02.0=

sin[sin

]

〖提示〗

●在适当范围内，作图观察一元复杂函数的形态：观察解的存在性；解的唯一性。进而，

借助图形法求近似解。

●匿名函数的使用方法。

●fzero指令的用法。

〖目的〗

●作图法求一元复杂函数解上的作用：观察解的存在性；解的唯一性；得近似解。

●匿名函数的使用方法。

●fzero指令的用法。

〖解答〗

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

matlab课后习题答案第四章

相关文档

最新文档