《第3章 数值计算-刘东卓》由会员分享，可在线阅读，更多相关《第3章 数值计算-刘东卓（74页珍藏版）》请在文档大全上搜索。
1、1矩阵的生成1.命令窗口直接输入法规则： 矩阵元素必须用 括住 矩阵元素必须用逗号或空格分隔 在 内矩阵的行与行之间必须用分号分隔2 矩阵元素可以是任何matlab表达式 ，可以是实数 ，也可以是复数，复数可用特殊函数i，j 输入 a=1 2 3;4 5 6 x=2 pi/2;sqrt(3) 3+5i 矩阵元素矩阵元素3符号的作用 逗号和分号的作用 分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。 分号如果出现在指令后，屏幕上将不显示结果。4注意：只要是赋过值的变量，不管是否在屏幕上显示过，都存储在工作空间中，以后可随时显示或调用。变量名尽可能不要重复，否则会覆盖 。 当一个指
2、令或矩阵太长时，可用续行5 冒号的作用 用于生成等间隔的向量，默认间隔为1。 用于选出矩阵指定行、列及元素。 循环语句62.用matlab函数创建矩阵 空阵 matlab允许输入空阵，当一项操作无结果时，返回空阵。 rand 随机矩阵 eye 单位矩阵 zeros 全部元素都为0的矩阵 ones 全部元素都为1的矩阵7 还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。 matlab函数名必须小写。83.由m文件生成94.由文本文件生成1031 矩阵的运算 矩阵的转置：用标号“ ”即可A = mag
3、ic(4) % 创建一个创建一个4阶魔方矩阵阶魔方矩阵B = A % 计算矩阵的转置计算矩阵的转置A = 16 2 3 13 5 11 10 8 9 7 6 12 4 14 15 1 B =16 5 9 4 2 11 7 14 3 10 6 15 13 8 12 111矩阵的对角元素 利用函数利用函数diag()进行对角运算进行对角运算A = pascal(3)A = 1 1 1 1 2 3 1 3 6v = diag(A)v = 1 2 6B = diag(v,1)B = 0 1 0 0 0 0 2 0 0 0 0 6 0 0 0 012矩阵的基本数学运算1. 矩阵加、减（,）运算规则： 相
4、加、减的两矩阵必须有相同的行和列两矩阵对应元素相加减。 允许参与运算的两矩阵之一是标量。标量与矩阵的所有元素分别进行加减操作。132. 矩阵乘（）运算规则： A矩阵的列数必须等于B矩阵的行数 标量可与任何矩阵相乘。a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*bc =14 32 23 14 3 矩阵除的运算在线性代数中没有，有矩阵逆的运算，在matlab中有两种矩阵除运算15稀疏矩阵 稀疏矩阵的生成 在MATLAB 7中，生成稀疏矩阵用特殊的函数来进行，这些函数有speye、spones、spdiags、sparse、find、full、spalloc、sprand和spra
5、ndn 等。16Speye & Sparse函数应用举例 A=eye(5)A = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 speye(size(A)ans = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (5,5) 1 sparse(A)ans = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (5,5) 117稀疏矩阵和全元素矩阵的转换 full(sparse(A)ans = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 sparse
6、(A)ans = (1,1) 1 (2,2) 1 (3,3) 1 (4,4) 1 (5,5) 1 A=eye(5)A = 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0 0 0 118稀疏带状矩阵 A = 11 0 13 0 0 22 0 24 0 0 33 0 41 0 0 44 0 52 0 0 0 0 63 0 0 0 0 74; B,d = spdiags(A)B = 41 11 0 52 22 0 63 33 13 74 44 24d = -3 0 2 B1 = magic(4)B1 = 16 2 3 13 5 11 10 8 9 7 6 12
7、 4 14 15 1 d1 = -5,-2,0,2; A = spdiags(B1,d1,A); full(A)ans = 3 0 12 0 0 10 0 1 2 0 6 0 41 11 0 15 0 52 7 0 16 0 63 14 0 5 0 7419说明说明 对稀疏矩阵进行操作，主要由对稀疏矩阵进行操作，主要由nnz、nonzeros、nzmax、sponse、spalloc、isspase、spyfun和和spy等函数来实现等函数来实现 。20 三角分解(lu) 正交分解(qr) 特征值分解(eig) Chollesky分解(chol) 奇异值分解(svd) 3.2 矩阵的分解21
8、行列式和逆矩阵 A = magic(3)A = 8 1 6 3 5 7 4 9 2%产生产生3阶魔方矩阵阶魔方矩阵 d = det(A)d = -360% 返回行列式的值返回行列式的值 B = inv(A)B = 0.1472 -0.1444 0.0639 -0.0611 0.0222 0.1056 -0.0194 0.1889 -0.1028%返回返回A的逆矩阵的逆矩阵22特征值分解 A = magic(3)A = 8 1 6 3 5 7 4 9 2 d = eig(A)d = 15.0000 4.8990 -4.8990%d为为A的特征值的特征值 V,D = eig(A)V = -0.57
9、74 -0.8131 -0.3416 -0.5774 0.4714 -0.4714 -0.5774 0.3416 0.8131D = 5.0000 0 0 0 4.8990 0 0 0 -4.8990% V为为A的广义特征值向量矩阵的广义特征值向量矩阵 D为为A的广义特征值矩阵的广义特征值矩阵23奇异值分解奇异值分解 A = 1 2;3 4;5 6;7 8A = 1 2 3 4 5 6 7 8 s = svd(A)s = 14.2691 0.6268 s=svd(X) U,S,V=svd(X) U,S,V=svd(X,0) %s为奇异值，为奇异值，X=USV24 U,S,V = svd(A,0
10、)U = -0.1525 -0.8226 -0.3499 -0.4214 -0.5474 -0.0201 -0.7448 0.3812S = 14.2691 0 0 0.6268V = -0.6414 0.7672 -0.7672 -0.6414 U,S,V = svd(A)U = 0.1525 0.8226 0.3945 0.3800 0.3499 0.4214 0.2428 0.8007 0.5474 0.0201 0.6979 0.4614 0.7448 0.3812 0.5462 0.0407S = 14.2691 0 0 0.6268 0 0 0 0V = -0.6414 0.767
11、2 -0.7672 -0.6414253.3 求解线性求解线性代数方程组方程组 x=Ab x=inv(A)*bA = 1 2 3;1 3 5;1 3 6;b = 2 4 5; % 将线性方程组系数改写成矩阵形式将线性方程组系数改写成矩阵形式x1 = Ab % 左除法求解方程组左除法求解方程组x2 = inv(A)*bx1 = -1 0 1x2 = -1 0 126Cholosky分解 返回正定矩阵的上三角矩阵返回正定矩阵的上三角矩阵A = pascal(3);b = 4 2 3;% 创建方程组矩阵，变量系数为帕斯卡矩阵创建方程组矩阵，变量系数为帕斯卡矩阵L = chol(A) % 计算计算A的
12、上三角矩阵的上三角矩阵x = L(Lb) % 求解方程组求解方程组A = 1 1 1 1 2 3 1 3 6L = 1 1 1 0 1 2 0 0 1x = 9 -8 327LU分解分解A = 2 1 1;2 3 2;2 3 4;b = 4 7 8; % 创建线性方程组矩阵创建线性方程组矩阵L,U = lu(A) % 计算计算L和和U矩阵矩阵x = U(Lb) % 求解方程组求解方程组L = 1 0 0 1 1 0 1 1 1U = 2 1 1 0 2 1 0 0 2x = 1.1250 1.2500 0.500028QR分解分解A = 1 2 3;4 5 6;7 8 9;10 11 12;b
13、 = 1 3 5 7; % 创建线性方程组矩阵Q,R = qr(A) % 计算Q和R矩阵x = R(Qb) % 求解方程组Q = -0.0776 -0.8331 0.5456 -0.0478 -0.3105 -0.4512 -0.6919 0.4704 -0.5433 -0.0694 -0.2531 -0.7975 -0.7762 0.3124 0.3994 0.3748R = -12.8841 -14.5916 -16.2992 0 -1.0413 -2.0826 0 0 -0.0000 0 0 0 x = 0.5000 0 0.1667293.4 求解非线性方程求解非线性方程 f(x)=0
14、options)x0,fsolve(funxx0),fsolve(funxfunction F=myfun(x)F=2*x(1)-x(2)-exp(-x(1); -x(1)+2*x(2)exp(-x(2);% 同一工作目录下，在同一工作目录下，在MATLAB命令窗口运行下列指令命令窗口运行下列指令x0 = -6;-5; % 给变量初始值给变量初始值x = fsolve(myfun,x0) % 求解方程，求解方程，options缺省，等缺省，等价形式还可以写成价形式还可以写成x = fsolve(myfun,x0)x = 0.5671 0.5671303.4 函数的零点函数的零点 一元函数的零点
15、一元函数的零点格式：格式：x=fzero(fun,x0) x=fzero(fun,x0,options) x,fval=fsolve() 一般通过作图观察，然后通过计算不断一般通过作图观察，然后通过计算不断地逼近地逼近31x = -5:0.1:5;f = x.3-2*x-5;xlabel(x);ylabel(f(x);plot(x,f,x,0) % 计算计算2附近的零点和该点的函数值附近的零点和该点的函数值xout,f = fzero(zero1,2)xout = 2.0946f = -8.8818e-01632多元函数的零点多元函数的零点例3-17% 创建二元函数x = -1:.01:1;y
16、 = x;X,Y = meshgrid(x,y); % 产生平面网点坐标Z = sin(X.2-Y);v = -0.2,0,0.2; % 选取三条个等高值，更可靠地判断0等高线的位置contour (X,Y,Z,v) % 绘制函数等高线图，如图3-2hold on;plot(x,0.4,0.6,y)hold off;xout = fsolve(zero2,0.6 0.4) % 计算零点值，x=xout(1)，y=xout(2)。fsolve的使用方法见3.4节f = sin(xout(1).2-xout(2)% 计算 的值 yout = fsolve(zero2,-0.5 -0.5)333.6
17、 函数的极值点函数的极值点 一元函数的极值点格式： x=fminbnd(fun,x1,x2) x=fminbnd(fun,x1,x2,options) x,fval=fminbnd()34多元函数的极值点无约束条件格式： x=fminsearch(fun,x0) x=fminsearch (fun,x0, options) x,fval=fminunc()有约束条件格式： x=fmincon(fun,x0,A,b) x,fval=fmincon()353.7 数值微积分 差分差分例例3-22A = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18;% 生
18、成生成1维矩阵维矩阵A1 = reshape(A,6,3) % 转换为转换为36维矩阵维矩阵B1 = diff(A1) % 求求1维维1阶差分阶差分B2 = diff(A1,1,2) % 求求2维维1阶差分阶差分B3 = diff(A1,2) % 求求1维维2阶差分阶差分A2 = reshape(A,2,3,3) % 转换为转换为233维矩阵维矩阵A2(:,:,1) =B4 = diff(A2,1,3) % 求求3维维1阶差分阶差分36梯度和偏导数格式：FX=gradient(F) FX,FY=gradient(F) =gradient(F,h)功能：求得函数F的梯度，FX,FY表示沿x，y方
19、向的偏导数37例3-23% 生成二元函数v = -2:0.2:2;x,y = meshgrid(v);z = x .* exp(-x.2 - y.2); % 绘制曲面，如图绘制曲面，如图3-4所示所示figure(1)mesh(x,y,z);px,py = gradient(z,.2,.2); % 求偏导数求偏导数figure(2)contour(v,v,z) % 绘制等高线，如图绘制等高线，如图3-5所示所示hold onquiver(v,v,px,py) % 绘制矢量场图，小箭头表示梯度绘制矢量场图，小箭头表示梯度hold off38一元函数的数值积分问题Sum(2*y(1:end-1,:
20、)+diff(y).*diff(x)/2由给定数据进行梯形求积由给定数据进行梯形求积39 格式： S=trapz(x,y) 例： x1=0:pi/30:pi; y=sin(x1) cos(x1) sin(x1/2); x=x1 x1 x1; S=sum(2*y(1:end-1,:)+diff(y).*diff(x)/2S = 1.9982 0.0000 1.9995 S1=trapz(x1,y) % 得出和上述完全一致的结果S1 = 1.9982 0.0000 1.999540单变量数值积分问题求解 格式： y=quad(Fun,a,b) y=quadl(Fun,a,b) % 求定积分 y=q
21、uad(Fun,a,b, ) y=quadl(Fun,a,b, ) %限定精度的定积分求解，默认精度为106。41 例：第三种：匿名函数(MATLAB 7.0)第二种：inline 函数第一种，一般函数方法42函数定义被积函数： y=quad(c3ffun,0,1.5)y = 0.9661 用 inline 函数定义被积函数： f=inline(2/sqrt(pi)*exp(-x.2),x); y=quad(f,0,1.5)y = 0.9661 运用符号工具箱： syms x, y0=vpa(int(2/sqrt(pi)*exp(-x2),0,1.5),60) y0 = .9661051464
22、75310713936933729949905794996224943257461473285749 y=quad(f,0,1.5,1e-20) % 设置高精度，但该方法失效43 例：提高求解精度： y=quadl(f,0,1.5,1e-20)y = 0.9661 abs(y-y0)ans = .6402522848913892e-16 format long 16位精度 y=quadl(f,0,1.5,1e-20)y = 0.9661051464753144 例：求解绘制函数: x=0:0.01:2, 2+eps:0.01:4,4; y=exp(x.2).*(x2); y(end)=0; x
23、=eps, x; y=0,y; fill(x,y,g)45 调用quad( ): f=inline(exp(x.2).*(x2)./(4-sin(16*pi*x),x); I1=quad(f,0,4)I1 = 57.76435412500863 调用quadl( ): I2=quadl(f,0,4)I2 = 57.76445016946768 syms x; I=vpa(int(exp(x2),0,2)+int(80/(4-sin(16*pi*x),2,4) I = 57.764450125053010333315235385182464.1 微积分问题的解析解 4.1.1 极限问题的解析解
24、单变量函数的极限 格式1： L= limit( fun, x, x0) 格式2： L= limit( fun, x, x0, left 或 right)47多重数值积分 二重积分函数如下：q=dblquad(fun,xmin,xmax,ymin,ymax)q=dblquad(fun,xmin,xmax,ymin,ymax,tol)q=dblquad(fun,xmin,xmax,ymin,ymax,tol,method) 或用符号计算可得到高精度解或用符号计算可得到高精度解48 例： syms x y z int(int(int(4*x*z*exp(-x2*y-z2),x,0,1),y,0,pi
25、),z,0,pi)ans =(Ei(1,4*pi)+log(pi)+eulergamma+2*log(2)*pi2*hypergeom(1,2,-pi2)Ei(n,z)为指数积分，无解析解，但可求其数值解： vpa(ans,60) ans = 3.10807940208541272283461464767138521019142306317021863483588493.8 多项式的运算 多项式的根 n次多项式具有次多项式具有n个根，当然这些根可能个根，当然这些根可能是实根，也可能含有若干对共轭复根。是实根，也可能含有若干对共轭复根。MATLAB提供的提供的roots函数用于求多项式函数用于求
26、多项式的全部根，其调用格式为：的全部根，其调用格式为： x=roots(P) 其中其中P为多项式的系数向量，求得的根赋为多项式的系数向量，求得的根赋给向量给向量x，即，即x(1),x(2),x(n)分别代表多分别代表多项式的项式的n个根。个根。50例例 求多项式求多项式x4+8x3-10的根。的根。命令如下：命令如下：A=1,8,0,0,-10;x=roots(A)若已知多项式的全部根，则可以用若已知多项式的全部根，则可以用poly函函数建立起该多项式，其调用格式为：数建立起该多项式，其调用格式为：P=poly(x)若若x为具有为具有n个元素的向量，则个元素的向量，则poly(x)建立建立以以
27、x为其根的多项式，且将该多项式的系为其根的多项式，且将该多项式的系数赋给向量数赋给向量P。51多项式的求值多项式的求值 MATLAB提供了两种求多项式值的函数：提供了两种求多项式值的函数：polyval与与polyvalm，它们的输入参数均，它们的输入参数均为多项式系数向量为多项式系数向量P和自变量和自变量x。两者的。两者的区别在于前者是代数多项式求值，而后区别在于前者是代数多项式求值，而后者是矩阵多项式求值。者是矩阵多项式求值。521代数多项式求值代数多项式求值polyval函数用来求代数多项式的值，其调函数用来求代数多项式的值，其调用格式为：用格式为：Y=polyval(P,x)若若x为一
28、数值，则求多项式在该点的值；若为一数值，则求多项式在该点的值；若x为向量或矩阵，则对向量或矩阵中的每为向量或矩阵，则对向量或矩阵中的每个元素求其多项式的值。个元素求其多项式的值。例例 已知多项式已知多项式x4+8x3-10，分别取，分别取x=1.2和和一个一个23矩阵为自变量计算该多项式的矩阵为自变量计算该多项式的值。值。532矩阵多项式求值矩阵多项式求值polyvalm函数用来求矩阵多项式的值，其调用格函数用来求矩阵多项式的值，其调用格式与式与polyval相同，但含义不同。相同，但含义不同。polyvalm函数函数要求要求x为方阵，它以方阵为自变量求多项式的为方阵，它以方阵为自变量求多项式
29、的值。设值。设A为方阵，为方阵，P代表多项式代表多项式x3-5x2+8，那么，那么polyvalm(P,A)的含义是：的含义是：A*A*A-5*A*A+8*eye(size(A)而而polyval(P,A)的含义是：的含义是：A.*A.*A-5*A.*A+8*ones(size(A)例例 仍以多项式仍以多项式x4+8x3-10为例，取一个为例，取一个22矩阵矩阵为自变量分别用为自变量分别用polyval和和polyvalm计算该多项计算该多项式的值。式的值。54多项式的四则运算多项式的四则运算1多项式的加减运算多项式的加减运算2多项式乘法运算多项式乘法运算函数函数conv(P1,P2)用于求多
30、项式用于求多项式P1和和P2的乘积。的乘积。这里，这里，P1、P2是两个多项式系数向量。是两个多项式系数向量。例例 求多项式求多项式x4+8x3-10与多项式与多项式2x2-x+3的乘积。的乘积。a = 1 8 0 0 -10;b = 2 -1 3; % 生成多项式向量生成多项式向量c = conv(a,b) % 计算乘积计算乘积s = poly2str(c,x) % 标准形式表示标准形式表示553多项式除法多项式除法函数函数Q,r=deconv(P1,P2)用于对多项式用于对多项式P1和和P2作除法运算。其中作除法运算。其中Q返回多项式返回多项式P1除以除以P2的商式，的商式，r返回返回P1
31、除以除以P2的余式。的余式。这里，这里，Q和和r仍是多项式系数向量。仍是多项式系数向量。deconv是是conv的逆函数，即有的逆函数，即有P1=conv(P2,Q)+r。56例例 求多项式求多项式x4+8x3-10除以多项式除以多项式2x2-x+3的结果的结果v = 1 8 0 0 -10;u = 2 1 3; % 生成多项式向量生成多项式向量c = conv(v,u) % 计算多项式乘积计算多项式乘积q1,r1 = deconv(c,v) % 求商，整除求商，整除r1为为0，商多，商多项式与项式与u相同相同q2,r2 = deconv(v,u) % 多项式求商，带余数多项式求商，带余数57
32、多项式部分分式展开b = 5 3 -2 7;a = -4 0 8 3; % 生成分子分母向量生成分子分母向量r,p,k = residue(b,a) % 部分分式展开部分分式展开b2,a2 = residue(r,p,k) % 进行部分分式逆运算进行部分分式逆运算r = -1.4167 -0.6653 1.3320p = 1.5737 -1.1644 -0.4093b2 = -1.2500 -0.7500 0.5000 -1.7500a2 = 1.0000 -0.0000 -2.0000 -0.7500k = -1.250058多项式的拟合-最小二乘拟合 格式： p=polyfit(x,y,n
33、) 功能：把自变量x和函数值y拟合成n阶多项式x0=0:0.1:1；y0=-.447 1.978 3.11 5.25 5.02 4.66 4.01 4.58 3.45 5.35 9.22;p=polyfit(x0,y0,3)p = 56.6915 -87.1174 40.0070 -0.9043xx=0:0.01:1;yy=polyval(p,xx);plot(xx,yy,-b,x0,y0,or)5960多项式的插值插值插值 插值的定义插值的定义是对某些集合给定的数据点之是对某些集合给定的数据点之间函数的估值方法。间函数的估值方法。 当不能很快地求出所需中间点的函数时，插值当不能很快地求出所需
34、中间点的函数时，插值是一个非常有价值的工具。是一个非常有价值的工具。 Matlab提供了一维、二维、提供了一维、二维、 三次样条等许多三次样条等许多插值选择插值选择61interp1 interp2 interp3 interpftinterpnspline griddata v 利用已知点确定未知点v 粗糙 精确v 集合大的 简化的插值函数插值函数62一维插值一维插值x = 0:10; % 给出已知基准数据给出已知基准数据y = sin(x);xi = 0:.25:10; % 在两个基准数据点间插入在两个基准数据点间插入3个点个点yi1 = interp1(x,y,xi,*nearest);
35、 yi2 = interp1(x,y,xi,*linear);yi3 = interp1(x,y,xi,*spline);yi4 = interp1(x,y,xi,*cubic); % 分别用四种方法求中间分别用四种方法求中间插值插值plot(x,y,o,xi,yi1,:,xi,yi2,-,xi,yi3,k.-,xi,yi4,-)% 绘图，并标注各曲线代表的插值方法，如图所示绘图，并标注各曲线代表的插值方法，如图所示legend(原始数据原始数据,最近点插值最近点插值,线性插值线性插值,样条插值样条插值,立方插值立方插值)6364二维插值二维插值axis(-3 3 -3 3 -5 20) %
36、确定坐标轴的范围确定坐标轴的范围X,Y = meshgrid(-3:.5:3); % 生成网格生成网格Z = peaks(X,Y); % 计算计算peaks函数值函数值figure(1)mesh(X,Y,Z) % 绘制原始数据曲面图，如图绘制原始数据曲面图，如图3-9 XI,YI = meshgrid(-3:.125:3); % 确定插值点确定插值点ZI1 = interp2(X,Y,Z,XI,YI,*nearest);% 计算用计算用nearest方法所得二维插值方法所得二维插值figure(2)mesh(XI,YI,ZI1) % 绘制最近点插值法曲面图，如绘制最近点插值法曲面图，如图图3-
37、10 65ZI2 = interp2(X,Y,Z,XI,YI,*linear);% 计算用计算用linear法所得法所得二维插值二维插值figure(3)mesh(XI,YI,ZI2) % 绘制双线性插值法曲面图，如图绘制双线性插值法曲面图，如图3-11ZI3 = interp2(X,Y,Z,XI,YI,*spline); % 计算用计算用spline法所得二维插值法所得二维插值figure(4)mesh(XI,YI,ZI3) % 绘制双样条插值法曲面图，如图绘制双样条插值法曲面图，如图3-12ZI4 = interp2(X,Y,Z,XI,YI,*cubic); % 计算用计算用cubic法所
38、得二维插值法所得二维插值figure(5)mesh(XI,YI,ZI4) % 绘制双立方插值法曲面图，如图绘制双立方插值法曲面图，如图3-136667686970713.9 初值常微分方程的求解 函数如下： ode23 ode45 ode113 ode23t ode15s ode23s ode23tb 微分方程数值解的函数格式微分方程数值解的函数格式t,Y=solver(odefun,tspan,y0,options)72例3-39先编写函数rigid.mfunction dy=rigid(t,y)dy=zeros(3,1);dy(1)=y(2)*y(3);Dy(2)=-y(1)*y(3);Dy(3)=-0.51*y(1)*y(2);73tspan = 0 12;y0 = 0 1 1; % 在同一目录下，计算方程数值解。在同一目录下，计算方程数值解。输入时间区间和初始条件输入时间区间和初始条件t,Y = ode45(rigid,tspan,y0);% 采用采用ode45算法求解方程，算法求解方程，options为默为默认值认值plot(t,Y(:,1),-,t,Y(:,2),-.,t,Y(:,3),.)% 绘制计算结果并标注，如图绘制计算结果并标注，如图3-14legend(Y1,Y2,Y3)74
