[数学实验与matlab]程序.doc

18数学实验与matlab 数学实验与matlab程序我将程序按书中的顺序排列，这样便于读者利用。 本书程序均通过了调式。可直接拷贝到命令窗口运行或编制m文件运行。如出现问题，可能是中英文标点的缘故（排版有可能使用了中文标点），请将中文标点换为英文标点试试。 实验矩阵运算与matlab命令1.1 知识要点与背景：知识要点和背景：日常矩阵及其运算【 a=4 2 3;1 3 2;1 3 3;3 2 2, % 表1-1、表1-2的数据分别写成矩阵形式 b=35 20 60 45;10 15 50 40;20 12 45 20 】【 c=a*b %矩阵乘法，求各订单所对应的原材料和劳动力 。 】 【 whos % 查看matlab工作空间中变量及其规模 】 1.2实验与观察：矩阵和matlab语言1.2.1 向量的生成和运算 【 x=linspace(0,4*pi,100); %将0,4区间100等分，产生了一个100维向量 y=sin(x); %计算函数值，产生了一个与x同维的100维函数向量y y1=sin(x).2; %计算函数向量，注意元素群运算 y2=exp(-x).*sin(x); %以x为横坐标，y为纵坐标画函数的图用不同的线型将函数曲线绘制在一个图上 plot(x,y,-,x,y1,-,x,y2,.-) 】 1. 向量的创建 直接输入向量。 【x1=1 2 4,x2=1,2,1,x3=x1 】 冒号创建向量 。【 x1=3.4:6.7 x2=3.4:2:6.7 x3=2.6:-0.8:0 】 生成线性等分向量。 【 x=linspace(0,1,5) 】2. 向量的运算 【 y=sin(x) 】 【 y1=sin(x).2; y2=exp(-x).*sin(x); 】 1.2.2.矩阵创建和运算1.创建矩阵（1）数值矩阵的创建 直接输入法创建简单矩阵。【 a=1 2 3 4; 5 6 7 8; 9 10 11 12 】 【 b=-1.3,sqrt(3);(1+2)*4/5,sin(5);exp(2),6 】 （2）符号矩阵的创建 【 syms a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34a1=a11 a12 a13 a14 ;a21 a22 a23 a24; a31 a32 a33 a34,b1=b11 b12 b13 b14 ;b21 b22 b23 b24; b31 b32 b33 b34 】 2.矩阵的运算 【 c=a1+b1,d=a1-b1 】 【 syms c ca=c*a1 】 【 c=a1*b1 】 ? error using = sym/mtimes, inner matrix dimensions must agree. 【 a2=a1(:,1:3), b1 】【 g=a2*b1 】【 g11=a2(1,:)*b1(:,1) 】 【 a, a_trans=a 】 【 h=1 2 3 ; 2 1 0 ; 1 2 3 , k=1 2 3 ; 2 1 0 ; 2 3 1h_det=det(h), k_det=det(k),h_inv=inv(h),k_inv=k-1 】 【 a=3 0 1; 1 1 0;0 1 4; b=inv(a-2*eye(3)*a, b=(a-2*eye(3)a 】3.分块矩阵：矩阵的裁剪、分割、修改与抽取（1）【 a=1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1, vr=1,3;vc=1,3;a1=a(vr,vc) %取出a的1、3行和1、3列的交叉处元素构成新矩阵a1 】 将上面的矩阵a分为四块，并把它们赋值到矩阵b中，观察运行后的结果。【 a11=a(1:2,1:2),a12=a(1:2,3:5),a21=a(3:4,1:2),a22=a(3:4,3:5)b=a11 a12;a21 a22 】a = 2 0 5 4 2 1 0 -1 2b = 1 2 4 -1 5 3 1 0 -1 0 2 3c = -3 4 18 13 13 14 20 -1 -7 -3 3 6 （2）矩阵的修改和提取 【 a=1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1 a(1,:)=0 0 0 0 0; a 】 观察：【 b(:,2,4)= %删除矩阵b的第2、4列 】（3）矩阵元素的抽取4.生成特殊矩阵 。 【 y1=rand(1,5), y2=rand(1,5),rand(seed,3), x1=rand(1,5), rand(seed,3), x2=rand(1,5) 】5. 常用矩阵函数6. 数据的简单分析 【 rand(seed,1);a=rand(3,6), asort=sort(a), amax=max(a), asum=sum(a) 】1.2.3 matlab工作环境和编程2.matlab的基本设计1.3应用、思考与练习1.3.1 关系矩阵1.3.2 投入产出 1.3.3 循环比赛的名次 【 a=0 1 1 0; 0 0 1 1; 0 0 0 1; 1 0 0 0, e=ones(4,1); c=a*e; s=c 】 画矩阵结构图的gplot指令。 （3） 【 clf, a=0 1 1 0;0 0 1 1;0 0 0 1;1 0 0 0; xy=0 1;0 0;-1 0.5;1 0.5; graphy_plot(a,xy,1,0.5), % gplot(a,xy) 】1.3.4 参考程序graphy_plot.m 【 function y=graphy_plot(a,xy,l,p) %画矩阵的有向结构图。 a为邻接矩阵，xy为顶点坐标，l控制参数，l=0，画无向图； %l=0,画有向图。p为控制箭头大小的参数。a=-max(abs(xy(:,1)*1.1;b=max(abs(xy(:,1)*1.1;c=-max(abs(xy(:,2)*1.1;d=max(abs(xy(:,2)*1.1;if l=0 gplot(a,xy),axis(a b c d),hold on,elseif l=0 u=;v=;x=;y=; n=length(a(:,1) ; for i=1:n k=find(a(i,:)=0); m=length(k); if(m=0) for j=1:m u(1)=(xy(k(j),1)-xy(i,1); v(1)=(xy(k(j),2)-xy(i,2); u(2)=eps; v(2)=eps; u=u;u;v=v;v; x=xy(i,1) xy(k(j),1);x; y=xy(i,2) xy(k(j),2);y; end text(xy(i,1),xy(i,2),bulletleftarrowfontsize16itv, um2str(i); hold on, end end gplot(a,xy),axis(a b c d),hold on, h=quiver(x,y,u,v,p);set(h,color,red);hold on, plot(xy(:,1),xy(:,2),k.,markersize,12),hold on,end , hold off 】实验函数的可视化与matlab作2.1 实验与观察：函数的可视化2.1.1 matlab二维绘图命令1.周期函数与线性p周期函数 观察 ： 【 clf, x=linspace(0,8*pi,100);f=inline(sin(x+cos(x+sin(x);y1=sin(x+cos(x+sin(x); y2=0.2*x+sin(x+cos(x+sin(x);plot(x,y1,k:,x,y2,k-) legend(sin(x+cos(x+sin(x),0.2x+sin(x+cos(x+sin(x),2) 】 2. plot指令：绘制直角坐标的二维曲线3. 图形的属性设置和屏幕控制 【 h=plot(0:0.1:2*pi,sin(0:0.1:2*pi); grid on set(h,linewidth,5,color,red); set(gca,gridlinestyle,-,fontsize,16) 】 设置y坐标的刻度并加以说明，并改变字体的大小。【 h=plot(0:0.1:2*pi,sin(0:0.1:2*pi);grid on， set(gca,ytick,-1 -0.5 0 0.5 1), set(gca,yticklabel,a|b|c|d|e), set(gca,fontsize,20) 】4. 文字标注指令 【 plot(x,y1,b,x,y2,k-) , set(gca,fontsize,15,fontname,times new roman), %设置轴对象的字体为times % new roman，字体的大小为15 title( itperoid and linear peroid function); %加标题，注意文字用单引号 加上 %斜杠后可输入不同的设置，例如it表示花括号里的文字为斜体；如果有多项设置， %则可用连续输入。 xlabel(x from 0 to 8*pi itt); ylabel(ity); %说明坐标轴 text(x(49),y1(50)-0.4,fontsize15bulletleftarrowthe period function itf(x); %在坐标(x(49),y1(50)-0.4）处作文字说明， 各项设置用隔开。 %fontsize15bulletleftarrow的意义依次是：字体大小=15 画圆点 左箭头 text(x(14),y2(50)+1,fontsize15the linear period function itg(x)rightarrowbullet) %与上一语句类似，用右箭头 】图2.5 文字标注 观察指令legend和num2str的用法：在同一张图上画出， 这里, 并进行适当的标注。 zxy2_2.m 【 clf, t=0:0.1:3*pi;alpha=0:0.1:3*pi; plot(t,sin(t),r-);hold on; plot(alpha,3*exp(-0.5*alpha),k:); set(gca,fontsize,15,fontname,times new roman), xlabel(itt(deg);ylabel(itmagnitude);title( itsine wave and itae-alphaittwave); %注意alpha的意义text(6,sin(6),fontsize15the value itsin(t) at itt=6rightarrowbullet, horizontalalignment,right), %上面的语句是一整行，如果要写成两行，必须使用续行号 ，例如要在“ bullet,” %后换行，需写“bullet, ”后才能换行。 % horizontalalignment,right 表示箭头所指的曲线对象在 文字的右边。text(2,3*exp(-0.5*2),fontsize15bulletleftarrow the value of it3e-0.5 itt=,num2str(3*exp(-0.5*2), at itt =2 ); %num2str的用法：string1 ,num2str,string2，注意方括号的使用。%legend(itsin(t),itae-alphat) % 请结合图形观察此命令的使用 】 运行结果如图2.6所示。 5. 图形窗口的创建和分割 观察：【 clf，b=2*pi;x=linspace(0,b,50);for k =1:9 y=sin(k*x); subplot(3,3,k),plot(x,y),axis(0,2*pi,-1,1)end 】x=linspace(0,4*pi,150);y1=cos(x);y=cos(cos(cos(cos(cos(cos(cos(x)+sin(x)+sin(x)+sin(x);plot(x,y,x,y1,r-),grid 2.1.2多元函数的可视化与空间解析几何（三维图形） 本节通过高等数学的几个例子观察matlab的三维绘图功能和技巧。1. 绘制二元函数 观察：绘制 的图象，作定义域的裁剪。 （1）观察meshgrid指令的效果。【 a=-0.98;b=0.98;c=-1;d=1;n=10;x=linspace(a,b,n); y=linspace(c,d,n);x,y=meshgrid(x,y);plot(x,y,+) 】 三维绘图指令mesh、meshc、surf。 （2）做函数的定义域裁剪，观察上述三维绘图指令的效果。程序zxy2_4.m【 clear,clf,a=-1;b=1;c=-15;d=15;n=20;eps1=0.01;x=linspace(a,b,n);y=linspace(c,d,n);x,y=meshgrid(x,y);for i=1:n %计算函数值z ，并作定义域裁剪 for j=1:n if (1-x(i,j)eps1|x(i,j)-y(i,j)eps1 %if语句这样用 z(i,j)=nan; %作定义域裁剪，定义域以外的函数值为nan else z(i,j)=1000*sqrt(1-x(i,j)-1.*log(x(i,j)-y(i,j); end endend zz=-20*ones(1,n);plot3(x,x,zz),grid off,hold on %画定义域的边界线mesh(x,y,z) %绘图，读者可用meshz, surf, meshc在此替换之 xlabel(x),ylabel(y),zlabel(z), box on %把三维图形封闭在箱体里 】 运行了zxy2_4.m 以后，有关向量存储在工作空间中，在此基础上，观察上述等值线绘制指令的运行效果。 【 cs,h=contour(x,y,z,15); clabel(cs,h,labelspacing,244) 】2. 三元函数可视化: slice指令 观察： 绘制三元函数的可视化图形。【 clf,x=linspace(-2,2,40); y=x; z=x; x,y,z=meshgrid(x,y,z); w=x.2+y.2+z.2;slice(x,y,z,w,-1,0,1,-1,0,1,-1,0,1),colorbar 】3. 空间曲线及其运动方向的表现：plot3和quiver指令【 clf, t=0:0.1:1.5;vx=2*t;vy=2*t.2;vz=6*t.3-t.2;x=t.2;y=(2/3)*t.3;z=(6/4)*t.4-(1/3)*t.3; %由速度得到曲线plot3(x,y,z,r.-),hold on %画飞行轨迹 %算数值梯度，也就是重新计算数值速度矢量，这只是为了编程的方便，不是必须的图2.12 飞机的飞行轨迹与方向vx=gradient(x);vy=gradient(y);vz=gradient(z);quiver3(x,y,z,vx,vy,vz),grid on %画速度矢量图xlabel(x),ylabel(y),zlabel(z) 】2.2应用、思考和练习2.2.1 线性p周期函数zxy2_3_f.m【 function f=zxy2_3_f(x) f=sin(x+cos(x); 】zxy2_3.m【 clear，clfa=-8;b=12;n=300;xx=linspace(a,b,n);h=zxy2_3_f(xx);s(1)=0;for i=2:n s(i)=s(i-1)+quad(zxy2_3_f,xx(i-1),xx(i);endsubplot(1,2,1),plot(xx,s,k-),axis(a,b,-1.5,9)subplot(1,2,2),plot(xx,h;zeros(1,length(xx),k-),axis(a,b,-1.5,1.5) 】2.2.2 平面截割法和曲面交线的绘制用平行截面法讨论由曲面 构成的马鞍面形状。 zxy2_6.m （ 平行截割法示例 ， 本程序的绘制两曲面交线方法可以套用） 【 clf, a=-20;eps0=1;x,y=meshgrid(-10:0.2:10); %生成平面网格v=-10 10 -10 10 -100 100; %设定空间坐标系的范围colormap(gray) %将当前的颜色设置为灰色z1=(x.2-2*y.2)+eps; %计算马鞍面函数z1=z1(x,y)z2=a*ones(size(x); %计算平面 z2=z2(x,y)r0=abs(z1-z2)=eps0; %计算一个和z1同维的函数r0,当abs(z1-z2)eps0时，r0 =0。 %可用mesh(x,y,r0)语句观察它的图形，体会它的作用，该方法可以套用。zz=r0.*z2;xx=r0.*x;yy=r0.*y; %计算截割的双曲线及其对应的坐标subplot(2,2,2), %在第2图形窗口绘制双曲线 h1=plot3(xx(r0=0),yy(r0=0),zz(r0=0),+); set(h1,markersize,2),hold on,axis(v),grid onsubplot(2,2,1), %在第一图形窗口绘制马鞍面和平面 mesh(x,y,z1);grid,hold on;mesh(x,y,z2); h2=plot3(xx(r0=0),yy(r0=0),zz(r0=0),.); %画出二者的交线 set(h2,markersize,6),hold on,axis(v),for i=1:5 %以下程序和上面是类似的，通过循环绘制一系列的平面去截割马鞍面 a=70-i*30; %在这里改变截割平面 z2=a*ones(size(x); r0=abs(z1-z2)=1; zz=r0.*z2;yy=r0.*y;xx=r0.*x; subplot(2,2,3), mesh(x,y,z1);grid,hold on;mesh(x,y,z2);hidden off h2=plot3(xx(r0=0),yy(r0=0),zz(r0=0),.); axis(v),grid subplot(2,2,4), h4=plot3(xx(r0=0),yy(r0=0),zz(r0=0),o); set(h4,markersize,2),hold on,axis(v),grid onend 】 2.2.3 微分方程的斜率场 绘制微分方程 的斜率场，并将解曲线画在图中，观察斜率场和解曲线的关系。 zxy2_5.m ( 绘制一阶微分方程的斜率场和解曲线)【 clf,clear %清除当前所有图形窗口的图像，清除当前工作空间的内存变量。a=0;b=4;c=0;d=4;n=15;x,y=meshgrid(linspace(a,b,n),linspace(c,d,n); %生成区域中的网格。z=x.*y; %计算斜率函数。 fx=cos(atan(x.*y);fy=sqrt(1-fx.2); %计算切线斜率矢量。quiver(x,y,fx,fy,0.5),hold on,axis(a,b,c,d) %在每一网格点画出相应的斜率矢量，0.5是控制矢量大小的控制参数，可以调整。x,y=ode45(zxy2_5f,0,4,0.4); %求解微分方程。 %zxy2_5f.m是方程相应函数f(x,y)的程序，单独编制；x0,xs=0,4为求解区间； %y0=0.4为初始值；输出变量x,y分别为解轨线自变量和因变量数组。plot(x,y,r.-) %画解轨线 】 zxy2_5f.m （微分方程的函数子程序）【 function dy=zxy2_5f(x,y)dy=x.*y; 】 2.2.4颜色控制和渲染及特殊绘图指令1.地球表面的气温分布（sphere指令） 【 a,b,c=sphere(40);t=max(max(abs(c)-abs(c);surf(a,b,c,t);axis(equal),colormap(hot), shading flat,colorbar 】2.旋转曲面的生成：柱面指令cylinder和光照控制指令surfl 【 x=0:0.1:10;z=x;y=1./(x.3-2.*x+4); u,v,w=cylinder(y);surfl(u,v,w,45,45); shading interp 】3.若干特殊图形 运行下面程序，了解各指令的用法和效果。【 x=1:10; y=5 6 3 4 8 1 10 3 5 6;subplot(2,3,1),bar(x,y),axis(1 10 1 11)subplot(2,3,2),hist(y,x),axis(1 10 1 4)subplot(2,3,3),stem(x,y,k),axis(1 10 1 11)subplot(2,3,4),stairs(x,y,k), axis(1 10 1 11)subplot(2,3,5), x = 1 3 0.5 5;explode = 0 0 0 1;pie(x,explode)subplot(2,3,6),z=0:0.1:100; x=sin(z);y=cos(z).*10; comet3(x,y,z) 】实验函数式直接确定型模型3.1知识要点和背景：函数 直接确定性模型 3.2实验与观察：插值与拟合3.2.1 插值方法与多项式拟合的概念3.2.2 用matlab作插值和拟合3.2.3 用鼠标选节点 观察插值、拟合的效果 3.2.4 观察程序说明zxy3_1.m【 clf,a=-1;b=1;n=100; %用内联函数inline命令定义函数， %在后面可直接用于函数g的计算,要改变函数做实验，可按此格式重新定义gg=inline(x2-x4); xx=linspace(a,b,n);for i=1:n gx(i)=gxx(i); % 前面已经用inline命令定义了g，可以这样用g计算函数值endymin=min(gx)*0.8;ymax=max(gx)*1.2;%分四个界面画图g的图形，以便于结果比较subplot(2,2,1),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(近邻插值)subplot(2,2,2),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(线性插值)subplot(2,2,3),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(样条插值)subplot(2,2,4),plot(xx,gx,-),grid,hold on,axis(a b ymin ymax),title(多项式拟合) %用鼠标在屏幕上选点x,y,button = ginput(n)，可套用下面程序的格式button=1;x1=a;y1=gx(1);while button=1 xi,yi,button=ginput(1); subplot(2,2,1),h=plot(xi,yi,ro) %在4个图形窗口画点 subplot(2,2,2),h=plot(xi,yi,ro) subplot(2,2,3),h=plot(xi,yi,ro) subplot(2,2,4),h=plot(xi,yi,ro) x1=xi,x1;y1=yi,y1; %将选的点存于向量x1,y1end x1=b,x1;y1=gx(n),y1; xx=linspace(a,b,n); %定义自变量xx %计算不同的插值函数：x1,y1为节点，xx为输入自变量 ynearest=interp1(x1,y1,xx,nearest); ylinear=interp1(x1,y1,xx,linear); yspline=interp1(x1,y1,xx,spline); %多项式拟合指令p,s = polyfit(x,y,n)，n为拟合多项式次数，x,y为被拟合数据， %p为拟合多项式的系数，s是用来做误差 估计和预测的数据结构。 p,c=polyfit(x1,y1,4); ypolyfit=polyval(p,xx); %用polyval(p,x)计算系数为p的多项式在标量或向量x处的值 subplot(2,2,1),h=plot(xx,ynearest,r-);set(h,linewidth,2) %画图subplot(2,2,2),h=plot(xx,ylinear,r-);set(h,linewidth,2);subplot(2,2,3),h=plot(xx,yspline,r-);set(h,linewidth,2) subplot(2,2,4),h=plot(xx,ypolyfit,r-);set(h,linewidth,2) 】3.3应用、思考和练习3.3.1若干函数的插值和拟合练习3.3.2几个应用问题1. 机床加工和水深流速问题 2. 内燃机转角与升程的关系 3. 随高度变化的大气压强4. 交通事故的调查3.3.4多元函数的插值 plot(fi,x,o) fi=91 105 110 115 120 124 128; x=0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409; plot(fi,x,o) x = 0 1.0869 1.9710 2.7555 3.3986 4.9073 8.3409 1. 矩形温箱的温度2. 航行区域的警示线3.3.5 fourier级数和周期函数的经验公式 zxy3_2.m【 clf,clear,n=10;m=3;x=1:1:12;y=3.1 3.8 6.9 12.7 16.8 20.5 24.5 25.9 22.0 16.1 10.7 5.4;z=(pi/6)*x;plot(z(1:12),y(1:12),o);hold onk=1:m; %计算数据矩阵。for i=1:n x(i,:)=1 cos(k*z(i) sin(k*z(i);endc=inv(x*x)*x*y(1:n), %求解。z1=linspace(0,2*pi,30);s=; %开始计算f级数的部分和。for i=1:30; sd=1 cos(k*z1(i) sin(k*z1(i); %注意k是向量。 sd=c.*sd; s=s,sum(sd);endplot(z1,s,r-),hold on, xlabel(月份),ylabel(平均气温) 】3.3.6由实验到理论:从开普勒定律和牛顿万有引力定律x=0 3 5 7 9 11 12 13 14 15;y=0 1.2 1.7 2.0 2.1 2.0 1.8 1.2 1.0 1.6;x1=0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9;y1=3.1950, 3.2299, 3.2532, 3.2611, 3.2516, 3.2282, 3.1870, 3.1266, 3.0594, 2.9757;clf,subplot(1,2,1),plot(x,y,k.),axis(-0.5 16 -0.5 4),subplot(1,2,2),plot(x1,y1,k.) t=0.241 0.615 1.000 1.881 11.862 29.457; r=0.387 0.723 1.000 1.524 5.203 9.539; t=t.2, r=r.3 plot(t,r,o) 运行结果（图形略，显示为直线）t = 0.0581 0.3782 1.0000 3.5382 140.7070 867.7148r = 0.0580 0.3779 1.0000 3.5396 140.8515 867.977775实验. 微分、积分和微分方程4.1. 知识要点和背景：微积分学基本定理 对于定义在a,b区间上的函数y=f(x),将a,b区间n等分，步长h=(b-a)/n.设分点为a=x0x1x2x=1 1 1 ;1 1 1;1 1 1,i1=cumsum(x,1),i2=cumsum(x,2)x = 1 1 1 1 1 1 1 1 1i1 = 1 1 1 2 2 2 3 3 3i2 = 1 2 3 1 2 3 1 2 3 观察：用累积和命令cumsum计算积分。【 x=linspace(0,pi,50); y=sin(x);t=cumsum(y)*pi/(50-1);i=t(50) 】 观察梯形公式计算的效果。【 x=linspace(0,pi,50); y=sin(x);t=trapz(x,y) 】 观察辛普森积分公式的计算效果。【 i1=quad(sin,0,pi), i2=quad8(sin,0,pi), 】 观察：用解常微分方程的方法计算积分。【 y0=0;x,y=ode23(myfun2_2,0,pi,y0);s=length(y);y(s) 】【 y0=0;x,y=ode45(myfun2_2,0,pi,y0);s=length(y);y(s) 】4.2.3 观察程序及其说明 zxy4_1.m （观察方程的解曲线族，图4.1）【 clf,clear,a=0;b=pi;c=0;d=2.2;n=20;x,y=meshgrid(linspace(a,b,n),linspace(c,d,n);z=x.*y; fx=cos(atan(sin(x);fy=sin(atan(sin(x); quiver(x,y,fx,fy,0.5),hold on,axis(a,b,c,d)x,y=ode45(zxy4_1f,0,pi,0); x1,y1=ode45(zxy4_1f,0,pi,0.2);x2,y2=ode45(zxy4_1f,0,pi,0.4);x3,y3=ode45(zxy4_1f,0,pi,0.6); plot(x,y,r.-,x1,y1,r.-,x2,y2,r.-,x3,y3,r.-), xlabel(x) ,ylabel(y) 】zxy4_1f.m 【 function dy=zxy4_1f(x,y) dy=sin(x); 】. 观察 程序zxy4_2.m （图4.14.3）【 clf, a=0;b=4*pi;n=31; x=linspace(a,b,n); y=zeros(1,n); y(1)=0; s(1)=0; for i=1:n-1 dy=myfun2_2(x(i); y(i+1)=y(i)+dy*(x(i+1)-x(i); hh(i)=myfun2_2(x(i); s(i+1)=s(i)+hh(i)*(x(i+1)-x(i); enda1=0.9*a;b1=1.1*b; %