matlab课后习题答案3.pdf

《matlab课后习题答案3.pdf》由会员分享，可在线阅读，更多相关《matlab课后习题答案3.pdf（19页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1 1 数字 1.5e2，1.5e3 中的哪个与 1500 相同吗？1.5e32 2 请指出如下 5 个变量名中，哪些是合法的？abcd-2 xyz_3 3chan a 变量 ABCDefgh 2、5 是合法的。3 3 在 MATLAB环境中，比 1 大的最小数是多少？1+eps4 4 设 a=-8,运行以下三条指令，问运行结果相同吗？为什么？w1=a(2/3)w2=(a2)(1/3)w3=(a(1/3)2w1=-2.0000+3.4641i ;w2=4.0000 ;w3=-2.0000+3.4641i5 5 指令 clear,clf,clc各有什么用处？clear去除工作空间中所有的变量。c

2、lf去除当前图形。clc第二章第二章去除命令窗口中所有显示。1 1 说出以下四条指令产生的结果各属于哪种数据类型，是“双精度对象，还是“符号符号对象？3/7+0.1 双双;sym(3/7+0.1)符符;sym(3/7+0.1)符符;vpa(sym(3/7+0.1)符符;2 2 在不加专门指定的情况下，以下符号表达式中的哪一个变量被认为是自由符号变量.sym(sin(w*t),sym(a*exp(-X),sym(z*exp(j*th)symvar(sym(sin(w*t),1)w a zw a z3 3 1试写出求三阶方程x 44.5 0正实根的程序。注意：只要正实根，不要出现其他根。2试求二阶

3、方程x2axa2 0在a1reset(symengine)syms x positivesolve(x3-44.5)ans=(2(2/3)*89(1/3)/2223 0时的根。2求五阶方程x ax a 0的实根syms a positive%注意：关于 x 的假设没有去除solve(x2-a*x+a2)Warning:Explicit solution could not be found.In solve at 83ans=empty sym syms x clearsyms a positivesolve(x2-a*x+a2)ans=a/2+(3(1/2)*a*i)/2 a/2-(3(1/

4、2)*a*i)/24 4 观察一个数在此用记述在以下四条不同指令作用下的异同。a=,b=sym(),c=sym(,d ),d=sym()在此，分别代表具体数值 7/3,pi/3,pi*3(1/3)；而异同通过 vpa(abs(a-d),vpa(abs(b-d),vpa(abs(c-d)等来观察。理解准确符号数值的创立法。高精度误差的观察。1x=7/3x=7/3;a=x,b=sym(x),c=sym(x,d),d=sym(7/3),a=2.3333b=7/3c=12.3333333333333334813630699500209d=7/3v1=vpa(abs(a-d),v2=vpa(abs(b-

5、d),v3=vpa(abs(c-d)v1=0.0v2=0.0v3=0.000000000000000148029736616687566666666677887162x=pi/3x=pi/3;a=x,b=sym(x),c=sym(x,d),d=sym(pi/3),a=1.0472b=pi/3c=1.047197551196597631317786181171d=pi/3v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d)v1=0.0v2=0.0v3=0.000000000000000114836428279922167628066158185543x=

6、pi*3(1/3)x=pi*3(1/3);a=x,b=sym(x),c=sym(x,d),d=sym(pi*3(1/3)a=4.5310b=1275352044764433/281474976710656c=4.5309606547207899041040946030989d=pi*3(1/3)v1=vpa(abs(a-d),v2=vpa(abs(b-d),v3=vpa(abs(c-d)v1=0.00000000000000026601114166290944374842393221638v2=0.00000000000000026601114166290944374842393221638v

7、3=0.0000000000000002660111416629094726767991785515a115 求符号矩阵A a21a31a12a22a32a13a23的行列式值和逆，所得结果应采用“子表达式置换简洁化。a33理解 subexpr 指令。A=sym(a11 a12 a13;a21 a22 a23;a31 a32 a33)DA=det(A)IA=inv(A);IAs,d=subexpr(IA,d)A=a11,a12,a13 a21,a22,a23 a31,a32,a33DA=2a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31+a13*a

8、21*a32-a13*a22*a31IAs=d*(a22*a33-a23*a32),-d*(a12*a33-a13*a32),d*(a12*a23-a13*a22)-d*(a21*a33-a23*a31),d*(a11*a33-a13*a31),-d*(a11*a23-a13*a21)d*(a21*a32-a22*a31),-d*(a11*a32-a12*a31),d*(a11*a22-a12*a21)d=1/(a11*a22*a33-a11*a23*a32-a12*a21*a33+a12*a23*a31+a13*a21*a32-a13*a22*a31)1k6 6 求x的符号解，并进而用该符号

9、解求()3k0k0k，()k01kk3，k0的准确值。symsum,subs 的应用。kx从实例中，感受指令所给出的关于k0符号解的含义。syms x kf=x(k);Z1=symsum(f,k,0,inf)Z1=piecewise(1=x,Inf,abs(x)1,-1/(x-1)subs(Z1,x,sym(-1/3),sym(1/pi),sym(3)ans=3/4,-1/(1/pi-1),Inf2 x 17 7 对于x 0，求k02k 1x12k1。提示：理论结果为ln x符号变量的限定性定义的作用。syms k;x=sym(x,positive);f_k=2/(2*k+1)*(x-1)/(

10、x+1)(2*k+1);s=simple(symsum(f_k,k,0,inf)%结果与理论值 lnx 相符！s=piecewise(abs(x-1)x+1,log(x)注意解答中，条件 abs(x-1)x+1 意味着：约束一：x-10此式总成立，说明“无约束。情况二：-(x-1)0此为“约束，满足题意。8 8 1通过符号计算求y(t)sint的导数dydt。2然后根据此结果，求dydt和t 0dydtt2。diff,limit指令的应用。如何理解运行结果。syms ty=abs(sin(t)d=diff(y)%求 dy/dtd0_=limit(d,t,0,left)%求 dy/dt|t=0-

11、dpi_2=limit(d,t,pi/2)%求 dy/dt|t=pi/2y=abs(sin(t)d=3sign(sin(t)*cos(t)d0_=-1dpi_2=09 9 求出1.710e xsin x dx的具有 64 位有效数字的积分值。符号积分的解析解和符号数值解。符号计算和数值计算的相互校验。1符号积分syms x clearsyms xy=exp(-abs(x)*abs(sin(x)si=vpa(int(y,-10*pi,1.7*pi),64)y=abs(sin(x)/exp(abs(x)si=1.0878494994129049131666718759481745208954585

12、352128459875194141662数值计算复验xx=-10*pi:pi/100:1.7*pi;sn=trapz(exp(-abs(xx).*abs(sin(xx)*pi/100sn=1.08771010 计算二重积分 12x21(x2 y2)dydx。变上限二重积分的符号计算法。syms x yf=x2+y2;r=int(int(f,y,1,x2),x,1,2)r=1006/1051111 在0,2区间，画出y(x)x0sin tdtt曲线，并计算y(4.5)。在符号计算中，经常遇到计算结果是特殊经典函数的情况。如何应用 subs 获得超过 16 位有效数字的符号数值结果。初步尝试 e

13、zplot 指令的简便。1符号计算syms t x;f=sin(t)/t;y=int(f,t,0,x)%将得到一个特殊经典函数y5=subs(y,x,sym(4.5)ezplot(y,0,2*pi)y=sinint(x)y5=1.65414041437924398350392248685154sinint(x)1.81.61.41.210.80.60.40.200123x4562数值计算复验tt=0:0.001:4.5;tt(1)=eps;yn=trapz(sin(tt)./tt)*0.001yn=1.65411212 在n 0的限制下，求y(n)20sinnxdx的一般积分表达式，并计算y(

14、1)的 32 位有效3数字表达。一般符号解与高精度符号数值解。syms xsyms n positivef=sin(x)n;yn=int(f,x,0,pi/2)y3s=vpa(subs(yn,n,sym(1/3)y3d=vpa(subs(yn,n,1/3)yn=beta(1/2,n/2+1/2)/2y3s=1.2935547796148952674767575125656y3d=1.29355477961489517824134054535531313 求方程x2 y21,xy 2的解。solve指令中，被解方程的正确书写，输出量的正确次序。eq1=x2+y2=1;eq2=x*y=2;x,y=

15、solve(eq1,eq2,x,y)x=(1/2+(15(1/2)*i)/2)(1/2)/2-(1/2+(15(1/2)*i)/2)(3/2)/2-(1/2+(15(1/2)*i)/2)(1/2)/2+(1/2+(15(1/2)*i)/2)(3/2)/2 (1/2-(15(1/2)*i)/2)(1/2)/2-(1/2-(15(1/2)*i)/2)(3/2)/2-(1/2-(15(1/2)*i)/2)(1/2)/2+(1/2-(15(1/2)*i)/2)(3/2)/2y=(1/2+(15(1/2)*i)/2)(1/2)-(1/2+(15(1/2)*i)/2)(1/2)5 (1/2-(15(1/2

16、)*i)/2)(1/2)-(1/2-(15(1/2)*i)/2)(1/2)14 求微分方程yyx 0的通解，并绘制任意常数为 1 时解的图形。54理解指令 dsolve 的正确使用。对 dsolve 输出结果的正确理解。ezplot 指令绘图时，如何进行线色控制。如何覆盖那些不能反映图形窗内容的图名。1求通解reset(symengine)clearsyms y xy=dsolve(0.2*y*Dy+0.25*x=0,x)y=2(1/2)*(C3-(5*x2)/8)(1/2)-2(1/2)*(C3-(5*x2)/8)(1/2)2根据所得通解中不定常数的符号写出“对其进行数值替代的指令yy=su

17、bs(y,C3,1)%将通解中的 C3 用 1 代替yy=2(1/2)*(1-(5*x2)/8)(1/2)-2(1/2)*(1-(5*x2)/8)(1/2)3观察通解中两个分解的平方是否相同yy(1)2=yy(2)2ans=14于是可考虑函数的平方关系syms Yfxy=Y2-yy(1)2fxy=Y2+(5*x2)/4-25根据平方关系式画完整曲线clfezplot(fxy,-2,2,-2,2)axis squaregrid onY2+(5 x2)/4-2=021.510.50-0.5-1-1.5-2-2x-1.5-1-0.56假设直接用“分解画曲线，那么将是不完整的0Y0.511.526 e

18、zplot(yy(1),hold oncc=get(gca,Children);set(cc,Color,r)ezplot(yy(2),axis(-2 2-2 2)legend(y(1),y(2),hold off;title()%覆盖不完全的图名gridaxis square2y(1)1.5y(2)10.50-0.5-1-1.5-2-2-1.5-1-0.500.511.52x1515 求一阶微分方程x at2bt,x(0)2的解。初值微分方程的符号解。pretty指令的使用。x=dsolve(Dx=a*t2+b*t,x(0)=2,t)pretty(x)%比拟易读的表达形式x=(t2*(3*b

19、+2*a*t)/6+2 2 t (3 b+2 a t)-+2 61616 求边值问题dfdx 3f 4g,dgdx 4 f 3g,f(0)0,g(0)1的解。注意：相应的数值解法比拟复杂。边值微分方程的符号解。f,g=dsolve(Df=3*f+4*g,Dg=-4*f+3*g,f(0)=0,g(0)=1)f=sin(4*t)*exp(3*t)g=cos(4*t)*exp(3*t)（1）数值数组及其运算习题 3及解答76 6 要求在闭区间0,2上产生具有 10 个等距采样点的一维数组。试用两种不同的指令实现。第 1 章数值计算中产生自变量采样点的两个常用指令的异同。%方法一t1=linspace

20、(0,2*pi,10)%方法二t2=0:2*pi/9:2*pi%要注意采样间距的选择，如这里的 2*pi/9.t1=Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 Columns 8 through 10 4.8869 5.5851 6.2832t2=Columns 1 through 7 0 0.6981 1.3963 2.0944 2.7925 3.4907 4.1888 Columns 8 through 10 4.8869 5.5851 6.28321 1 由指令 rng(default),A=rand(3

21、,5)生成二维数组 A，试求该数组中所有大于 0.5 的元素的位置，分别求出它们的“全下标和“单下标。第 1 章数组下标的不同描述：全下标和单下标。第 1 章sub2ind,int2str,disp 的使用。第 1 章随机发生器的状态控制：保证随机数的可复现性。rng(default)A=rand(3,5)ri,cj=find(A0.5);id=sub2ind(size(A),ri,cj);ri=ri;cj=cj;disp()disp(大于 0.5 的元素的全下标)disp(行号 ,int2str(ri)disp(列号 ,int2str(cj)disp()disp(大于 0.5 的元素的单下标

22、)disp(id)A=0.8147 0.9134 0.2785 0.9649 0.9572 0.9058 0.6324 0.5469 0.1576 0.4854 0.1270 0.0975 0.9575 0.9706 0.8003大于 0.5 的元素的全下标行号 1 2 1 2 2 3 1 3 1 3列号 1 1 2 2 3 3 4 4 5 5大于 0.5 的元素的单下标 1 2 4 5 8 9 10 12 13 151 1 矩阵12，运行指令 B1=A.(0.5),B2=A(0.5),可以观察到不同运算方法所得结果不同。1A 34数组运算和矩阵运算的不同。如何判断两个双精度数组是否相等。请分

23、别写出根据 B1,B2 恢复原矩阵 A的程序。2用指令检验所得的两个恢复矩阵是否相等。第 1 章第 1 章8第 1 章norm指令的应用。A=1,2;3,4;B1=A.0.5B2=A0.5A1=B1.*B1;A2=B2*B2;norm(A1-A2,fro)%求误差矩阵的 F-范数，当接近 eps 量级时，就认为实际相等B1=1.0000 1.4142 1.7321 2.0000B2=0.5537+0.4644i 0.8070-0.2124i 1.2104-0.3186i 1.7641+0.1458ians=8.4961e-0161 1 在时间区间 0,10中，绘制y 1e0.5tcos2t曲线

24、。要求分别采取“标量循环运算法和“数组运算法编写两段程序绘图。第 1 章加强理解数组运算的机理和应用。第 1 章初步使用subplot,plot,xlabel,ylabel等指令绘图。%标量循环运算法t=linspace(0,10,200);N=length(t);y1=zeros(size(t);for k=1:Ny1(k)=1-exp(-0.5*t(k)*cos(2*t(k);endsubplot(1,2,1),plot(t,y1),xlabel(t),ylabel(y1),grid on%数组运算法y2=1-exp(-0.5*t).*cos(2*t);subplot(1,2,2),plo

25、t(t,y2),xlabel(t),ylabel(y2),grid on1.51.511y10.5y20.5005t10005t101 1 先运行 clear,format long,rng(default),A=rand(3,3)，然后根据 A写出两个矩阵：一个对角阵 B，其相应元素由 A的对角元素构成；另一个矩阵 C，其对角元素全为 0，而其余元素与对应的 A阵元素相同。第 1 章clear,format longrng(default)A=rand(3,3)B=diag(diag(A)C=A-BA=常用指令 diag的使用场合。9 0.814723686393179 0.91337585

26、6139019 0.278498218867048 0.905791937075619 0.632359246225410 0.546881519204984 0.126986816293506 0.097540404999410 0.957506835434298B=0.814723686393179 0 0 0 0.632359246225410 0 0 0 0.957506835434298C=0 0.913375856139019 0.278498218867048 0.905791937075619 0 0.546881519204984 0.126986816293506 0.09

27、7540404999410 01 1先运行指令 x=-3*pi:pi/15:3*pi;y=x;X,Y=meshgrid(x,y);warning off;Z=sin(X).*sin(Y)./X./Y;产生矩阵Z。1请问矩阵 Z 中有多少个“非数数据？2用指令 surf(X,Y,Z);shading interp观察所绘的图形。3请写出绘制相应的“无裂缝图形的全部指令。第 1 章初步感受三维曲面的绘制方法。第 1 章非数 NaN的产生，非数的检测，和对图形的影响。第 1 章sum的应用。第 1 章eps 如何克服“被零除的为难。x=-3*pi:pi/15:3*pi;y=x;X,Y=meshgri

28、d(x,y);warning offZ=sin(X).*sin(Y)./X./Y;NumOfNaN=sum(sum(isnan(Z)%计算“非数数目subplot(1,2,1),surf(X,Y,Z),shading interp,title(有缝图)%产生无缝图XX=X+(X=0)*eps;YY=Y+(Y=0)*eps;ZZ=sin(XX).*sin(YY)./XX./YY;subplot(1,2,2),surf(XX,YY,ZZ),shading interp,title(无缝图)NumOfNaN=1811k 19k 12k 29k 2，当k依次1 1 下面有一段程序，企图用来解决如下计算

29、任务：有矩阵A Akk2k10k取 10,9,8,7,6,5,4,3,2,1 时，计算矩阵A Ak“各列元素的和，并把此求和结果存放为矩阵 Sa 的第 k行。例如k 3时，A1428阵为2529，此时它各列元素的和是一个(110)行数组36301061587，并把它保存为 Sa 的第 3 行。问题：该段程序的计算结果对吗？假设计算结果不正确，请指出错误发生的根源，并改正之。第 1 章正确理解 sum的工作机理。第 1 章reshape 的应用。1企图用以下程序完成题目要求。for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);endSaSa=55

30、 55 55 55 55 55 55 55 55 55 3 7 11 15 19 23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612

31、693 774 55 155 255 355 455 555 655 755 855 9552正确性分析除 k=1外，计算所得 Sa 所有行的结果都正确。但 k=1时，A A1 1,2,10，Sa的第 1 行应该与A A1相同。上述程序的错误是对 sum理解不正确。sum对二维数组，求和按列施行；而对一维数组，不管行数组或列数组，总是求那数组所有元素的和。正确的程序应该写成for k=10:-1:1A=reshape(1:10*k,k,10);Sa(k,:)=sum(A);if k=1Sa(k,:)=A;endendSaSa=1 2 3 4 5 6 7 8 9 10 3 7 11 15 19

32、23 27 31 35 39 6 15 24 33 42 51 60 69 78 87 10 26 42 58 74 90 106 122 138 154 15 40 65 90 115 140 165 190 215 240 21 57 93 129 165 201 237 273 309 345 28 77 126 175 224 273 322 371 420 469 36 100 164 228 292 356 420 484 548 612 45 126 207 288 369 450 531 612 693 774 55 155 255 355 455 555 655 755 85

33、5 955第 1 章数值运算习题 4 及解答 根据题给的模拟实际测量数据的一组t和然后把获得强调：要非常慎用数值导数计算。y(t)试用数值差分 diff或数值梯度 gradient 指令计算y(t)，prob_data401.maty(t)和y(t)曲线绘制在同一张图上，观察数值求导的后果。模拟数据从11练习 mat数据文件中数据的获取。实验数据求导的后果把两条曲线绘制在同一图上的一种方法。1从数据文件获得数据的指令假设 prob_data401.mat文件在当前目录或搜索路径上clearload prob_data401.mat2用 diff求导的指令dt=t(2)-t(1);yc=diff

34、(y)/dt;%注意 yc 的长度将比 y 短 1plot(t,y,b,t(2:end),yc,r)grid on1.510.50-0.5-1-1.5-2012345673用 gradent 求导的指令图形与上相似dt=t(2)-t(1);yc=gradient(y)/dt;plot(t,y,b,t,yc,r)grid on说明不到万不得已，不要进行数值求导。假假设一定要计算数值导数，自变量增量 dt 要取得比原有数据相对误差高 1、2 个量级以上。求导会使数据中原有的噪声放大。采用数值计算方法，画出y(x)x0sin tdtt在0,10区间曲线，并计算y(4.5)。提示指定区间内的积分函数可

35、用 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.6541y(4.5)1221.81.61.41.210.80.60.40.200246810 求函数f(x)esin3x的数值积分s f(x)dx，并请采用符号计算尝试复算。0提示数值积分均可尝试。符号积分的局限性。符

36、号积分的局限性。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)用 quad 求取1.7 5e xsin x dx的数值积分，并保证积分的绝对精度为109。quadl，精度可控，计算较快。近似积分指令 trapz 获得高精度积分的内存和时间代价较高。%精度可控的数值积分fx=(

37、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.0878494997343013%符号积分算法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

38、求函数f(t)(sin5t)2e0.06t1.5tcos2t 1.8t 0.5在区间5,5中的最小值点。2理解极值概念的邻域性。如何求最小值。学习运用作图法求极值或最小值。感受符号法的局限性。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？);%该指令只能在

39、指令窗中运行tt=linspace(-5,5,N+1);for k=1:Ntmin(k),fobj(k)=fminbnd(ft,tt(k),tt(k+1);endfobj,ii=sort(fobj);%将目标值由小到大排列tmin=tmin(ii);%使极小值点做与目标值相应的重新排列fobj,tmin2最后确定的最小值点在N 1,2,10的不同分割下，经观察，最后确定出最小值点是-1.28498111480531相应目标值是-0.186048010065453采用作图法近似确定最小值点另一方法A在指令窗中运行以下指令：clearft=(t)sin(5*t).2.*exp(0.06*t.*t)

40、+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,shgB经观察后，把最小值附近邻域放到足够大，然后运行以下指令，那放大图形被推向前台，与此同时光标变为“十字线，利用它点击极值点可得到最小值数据tmin2,fobj2=ginput(1)tmin2=-1.28500000993975fobj2=-0.1860479936913614出现具有相同数值的刻度区域说明已达最小可分辨状态4符号法求最小值的尝试syms tfts=sin(5*t)2*exp(0.06*t*t)-1.5*t*cos(2*t)+1.8

41、*abs(t+0.5);dfdt=diff(fts,t);%求导函数tmin=solve(dfdt,t)%求导函数的零点fobj3=subs(fts,t,tmin)%得到一个具体的极值点tmin=-.60100931947716486053884417850955e-2fobj3=.89909908144684551670208797723124说明最小值是对整个区间而言的，极小值是对邻域而言的。在一个区间中寻找最小值点，对不同子区间分割进行屡次搜索是必要的。这样可以防止把极小值点误作为最小值点。最小值点是从一系列极小值点和边界点的比拟中确定的。作图法求最小值点，很直观。假假设绘图时，自变量步

42、长取得足够小，那么所求得的最小值点有相当好的精度。符号法在本例中，只求出一个极值点。其余很多极值点无法秋初，更不可能得到最小值。d2y(t)dy(t)dy(0)3 2y(t)1,y(0)1,0，用数值法和符号法求y(t)t0.5。设2dtdtdt学习如何把高阶微分方程写成一阶微分方程组。ode45 解算器的导数函数如何采用匿名函数形式构成。如何从 ode45 一组数值解点，求指定自变量对应的函数值。1改写高阶微分方程为一阶微分方程组令y1(t)y(t),y2(t)dy(t)，于是据高阶微分方程可写出dtdy1(t)y2(t)dtdy(t)2 2y1(t)3y2(t)1dt2运行以下指令求 y(

43、t)的数值解format longts=0,1;y0=1;0;15dydt=(t,y)y(2);-2*y(1)+3*y(2)+1;%匿名函数写成的 ode45 所需得导数函数tt,yy=ode45(dydt,ts,y0);y_05=interp1(tt,yy(:,1),0.5,spline),%用一维插值求 y(0.5)y_05=0.789580207901273符号法求解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=.789580

44、35647060552916850705213780说明第条指令中的导数函数也可采用 M 函数文件表达，具体如下。function S=prob_DyDt(t,y)S=y(2);-2*y(1)+3*y(2)+1;求矩阵AxAx b b的解，A为 3 阶魔方阵，b 是(31)的全 1 列向量。提示了解 magic指令rref 用于方程求解。矩阵除法和逆阵法解方程。满秩方阵求解的一般过程。rref 用于方程求解。矩阵除法和逆阵法解方程。A=magic(3);%产生 3 阶魔方阵b=ones(3,1);%3*1全 1 列向量R,C=rref(A,b)%Gauss Jordan 消去法解方程，同时判断

45、解的唯一性x=Ab%矩阵除解方程xx=inv(A)*b%逆阵法解方程R=1.0000 0 0 0.0667 0 1.0000 0 0.0667 0 0 1.0000 0.0667C=1 2 3x=0.0667 0.0667 0.0667xx=0.0667 0.0667 0.0667说明rref指令通过对增广矩阵进行消去法操作完成解方程。由分解得到的 3 根“坐标向量和或C3 指示的 3 根基向量，可见 A3满秩，因此方程解唯一。在本例情况下，矩阵除、逆阵法、rref法所得解一致。求矩阵AxAx提示 b b的解，A为 4 阶魔方阵，b 是(41)的全 1 列向量。16用 rref 可观察 A不满

46、秩，但 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.000

47、0 -3.0000 -0.0588 0 0 0 0 0C=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=AbWarning

48、:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=1.306145e-017.x=0.0588 0.1176 -0.0588 0运行结果指示：矩阵除法给出的解与rref解相同。实际上，MATLAB在设计“除法程序时，针对“b 在 A值空间中的情况，就是用 rref求解的。3逆阵法的解xx=inv(A)*bWarning:Matrix is close to singular or badly scaled.Results may be inaccurate.RCOND=1.306145e-0

49、17.xx=0.0469 0.1875 -0.0625 -0.01563验证前面所得的解b_rref=A(:,C)*R(C,5)%验算 rref 的解b_d=A*x%验算矩阵除的解17b_inv=A*xx%验算逆阵法的解b_rref=1 1 1 1b_d=1 1 1 1b_inv=0.7344 1.5469 1.1719 1.8594显然，在本例中，逆阵法的解是错误的。原因是：A不满秩，A的逆阵在理论上不存在。这里所给出的逆阵是不可信的。4求齐次解xg=null(A)%Ax=0 的齐次解xg=0.2236 0.6708 -0.6708 -0.22365方程的全解齐次解的任何倍与特解之和就构成方

50、程的全解。下面通过一组随机系数验证。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.3077xx=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.27