第4章数值运算习题 4 及解答1 根据题给的模拟实际测量数据的一组t和)(t y试用数值差分diff或数值梯度gradient指令计算)(ty'曲线绘制y'，然后把)(t y和)(t 在同一张图上，观察数值求导的后果。
（模拟数据从prob_data401.mat获得）〖目的〗●强调：要非常慎用数值导数计算。
●练习mat数据文件中数据的获取。
●实验数据求导的后果●把两条曲线绘制在同一图上的一种方法。
〖解答〗（1）从数据文件获得数据的指令假如prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat（2）用diff求导的指令dt=t(2)-t(1);yc=diff(y)/dt; %注意yc的长度将比y短1plot(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 tt 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 ex f 3sin )(=的数值积分⎰=π0 )(dx x f s ，并请采用符号计算尝试复算。
〖提示〗● 数值积分均可尝试。
● 符号积分的局限性。
〖目的〗● 符号积分的局限性。
〖解答〗 dx=pi/2000;x=0:dx:pi;s=trapz(exp(sin(x).^3))*dxs =5.1370符号复算的尝试syms xf=exp(sin(x)^3);ss=int(f,x,0,pi)Warning: Explicit integral could not be found.> In sym.int at 58ss =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 longsq=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)))*dxst =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.0878494994129115 求函数5.08.12cos 5.1)5(sin )(206.02++-=t t t e t t f t 在区间]5,5[-中的最小值点。
〖目的〗● 理解极值概念的邻域性。
● 如何求最小值。
● 学习运用作图法求极值或最小值。
● 感受符号法的局限性。
〖解答〗（1）采用fminbnd 找极小值点在指令窗中多次运行以下指令，观察在不同数目子区间分割下，进行的极小值搜索。
然后从一系列极小值点中，确定最小值点。
clearft=@(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 ）在指令窗中运行以下指令：clearft=@(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.28500000993975fobj2 =-0.18604799369136出现具有相同数值的刻度区域表明已达最小可分辨状态（4）符号法求最小值的尝试syms tfts=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-2fobj3 =.89909908144684551670208797723124〖说明〗●最小值是对整个区间而言的，极小值是对邻域而言的。
●在一个区间中寻找最小值点，对不同子区间分割进行多次搜索是必要的。
这样可以避免把极小值点误作为最小值点。
最小值点是从一系列极小值点和边界点的比较中确定的。
●作图法求最小值点，很直观。
假若绘图时，自变量步长取得足够小，那么所求得的最小值点有相当好的精度。
●符号法在本例中，只求出一个极值点。
其余很多极值点无法秋初，更不可能得到最小值。
6 设0)0(,1)0(,1)(2)(3)(22===+-dtdy y t y dt t dy dt t y d ，用数值法和符号法求5.0)(=t t y 。
〖目的〗● 学习如何把高阶微分方程写成一阶微分方程组。
● ode45解算器的导数函数如何采用匿名函数形式构成。
● 如何从ode45一组数值解点，求指定自变量对应的函数值。
〖解答〗（1）改写高阶微分方程为一阶微分方程组令dtt dy t y t y t y )()(),()(21==，于是据高阶微分方程可写出 ⎪⎩⎪⎨⎧++-==1)(3)(2)()()(21221t y t y dtt dy t y dt t dy（2）运行以下指令求y(t)的数值解format longts=[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 39 55 5417 47 4640 26 2732 34 3541 23 2249 15 148 58 59B2 =-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.3536B3 =0.3536 0.6270 0.39130.3536 -0.4815 -0.24580.3536 -0.3361 -0.10040.3536 0.1906 -0.04510.3536 0.0451 -0.19060.3536 0.1004 0.33610.3536 0.2458 0.48150.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 10 1 0 0 1 0 3 4 -3 -4 70 0 1 0 0 1 -3 -4 4 5 -70 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 00 0 0 0 0 0 0 0 0 0 0B2_A=rref([B2,A])B2_A =Columns 1 through 71.0000 0 0 -91.9239 -91.9239 -91.9239 -91.92390 1.0000 0 51.8459 -51.8459 -51.8459 51.84590 0 1.0000 9.8995 -7.0711 -4.2426 1.41420 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0Columns 8 through 11-91.9239 -91.9239 -91.9239 -91.923951.8459 -51.8459 -51.8459 51.8459-1.4142 4.2426 7.0711 -9.89950 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0B3_A=rref([B3,A])B3_A =Columns 1 through 71.0000 0 0 91.9239 91.9239 91.9239 91.92390 1.0000 0 42.3447 -38.1021 -33.8594 29.61680 0 1.0000 12.6462 -16.8889 -21.1315 25.37410 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 00 0 0 0 0 0 0Columns 8 through 1191.9239 91.9239 91.9239 91.923925.3741 -21.1315 -16.8889 12.646229.6168 -33.8594 -38.1021 42.34470 0 0 00 0 0 00 0 0 00 0 0 00 0 0 0〖说明〗●magic(n)产生魔方阵。