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

t x y x ⎰

=0sin )(在]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 =

3 求函数x e

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⎰--π

π的数值积分，并保证积分的绝对精度为910-。

〖目的〗

● 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