第4章数值运算
习题 4 及解答
11  根据题给的模拟实际测量数据的一组和 试用数值差分diff或数值梯度gradient指令计算，然后把和曲线绘制在同一张图上，观察数值求导的后果。（模拟数据从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')
grid on
 
（3）用gradent求导的指令（图形与上相似）
dt=t(2)-t(1);
yc=gradient(y)/dt;
plot(t,y,'b',t,yc,'r')
grid on 
matlab求导
〖说明〗
●不到万不得已，不要进行数值求导。
●假若一定要计算数值导数，自变量增量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==4.5);   
s45=s(ii)
s45 =
    1.6541
 
13  求函数的数值积分，并请采用符号计算尝试复算。
〖提示〗
●数值积分均可尝试。
●符号积分的局限性。
〖目的〗
●符号积分的局限性。
〖解答〗
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)
14  用quad求取的数值积分，并保证积分的绝对精度为。
〖目的〗
●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 =
.0878********
%符号积分算法
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
15  求函数在区间中的最小值点。
〖目的〗
●理解极值概念的邻域性。
●如何求最小值。
●学习运用作图法求极值或最小值。
●感受符号法的局限性。
〖解答〗
（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    ？');%该指令只能在指令窗中运行
    tt=linspace(-5,5,N+1);
    for k=1:N
    [tmin(k),fobj(k)]=fminbnd(ft,tt(k),tt(k+1));
    end
    [fobj,ii]=sort(fobj);    %将目标值由小到大排列
    tmin=tmin(ii);                %使极小值点做与目标值相应的重新排列
fobj,tmin
（2）最后确定的最小值点
        在的不同分割下，经观察，最后确定出
最小值点是              -1.28498111480531
相应目标值是        -0.186********545
（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.186********136 
出现具有相同数值的刻度区域表明已达最小可分辨状态
（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
〖说明〗
●最小值是对整个区间而言的，极小值是对邻域而言的。
●在一个区间中寻最小值点，对不同子区间分割进行多次搜索是必要的。这样可以避免把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比较中确定的。
●作图法求最小值点，很直观。假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。
●符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。
16  设，用数值法和符号法求。
〖目的〗
●学习如何把高阶微分方程写成一阶微分方程组。
●ode45解算器的导数函数如何采用匿名函数形式构成。
●如何从ode45一组数值解点，求指定自变量对应的函数值。
〖解答〗
（1）改写高阶微分方程为一阶微分方程组
    令，于是据高阶微分方程可写出
（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函数文件表达，具体如下。