1、数学实验与Matlab,http:/ 周 晓 阳 华中科技大学数学系 ,实 验 一,矩阵运算与Matlab命令,日常矩阵及其运算,矩阵应用实例: 榄球防护用品的生产管理,应用问题,一个工厂生产三种橄榄球用品 ： 防护帽、 垫肩、臀垫。 需要不同数量的：硬塑料 、 泡沫塑料 尼龙线 、 劳动力。 为监控生产，管理者对它们之间的关系十分关心。 为把握这些量的关系，他列出下面的表,原料产品关系表,订单,管理者接到四份订单如上表所示。 问应该如何计算每份订单所需的原材料，以便组织生产？,将表格写成矩阵形式,计 算,输入下面Matlab指令 A=4 2 3;1 3 2;1 3 3;3 2 2, B=35
2、 20 60 45;10 15 50 40;20 12 45 20 C=A*B 请自行计算观看结果,Matlab基本指令,向量的创建和运算,1. 直接输入向量,x1=1 2 4,x2=1,2,1,x3=x1 运行结果x1 = 1 2 4x2 = 1 2 1x3 =124,2.冒号创建向量,x1=3.4:6.7,x2=3.4:2:6.7x3=2.6:-0.8:0 运算结果 x1 =3.4000 4.4000 5.4000 6.4000 x2 =3.4000 5.4000 x3 =2.6000 1.8000 1.0000 0.2000,3.生成线性等分向量,指令x=linspace(a,b,n)
3、在a,b区间产生 n 个等分点(包括端点) x=linspace(0,1,5) 结果 x = 0 0.2500 0.5000 0.7500 1.0000,工作空间,在Matlab窗口创建向量后并运行后，向量就存在于工作空间，可以被调用。,向量的运算,设x=x1 x2 x3; y=y1 y2 y3;为两个三维向量，a,b为标量。 向量的数乘：a*x=a*x1 a*x2 a*x3 向量的平移： x+b=x1+b x2+b x3+b 向量和： x+y=x1+y1 x2+y2 x3+y3 向量差： x-y=x1-y1 x2-y2 x3-y3 数的乘幂： 如 a2,元素群运算(四则运算）,x.*y=x1
4、*y1 x2*y2 x3*y3 (元素群乘积) x./y=x1/y1 x2/y2 x3/y3 (元素群右除，右边的y做分母) x.y=y1/x1 y2/x2 y3/x3 (元素群左除，左边的x做分母) x.5=x15 x25 x35 (元素群乘幂) 2.x=2x1 2x2 2x3 (元素群乘幂) x.y=x1y1 x2y2 x3y3 (元素群乘幂),元素群运算（函数计算）,Matlab有许多内部函数，可直接作用于向量产生一个同维的函数向量。 x=linspace(0,4*pi,100);（产生100维向量x） y=sin(x); (y也自动为100维向量) y1=sin(x).2; y2=ex
5、p(-x).*sin(x); 观察结果,创建矩阵（数值矩阵的创建）,直接输入法创建简单矩阵。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 观察运行结果,创建矩阵（符号矩阵的创建）,用指令“syms”说明符号变量。 syms a11 a12 a13 a14 a21 a22 a23 a24 a31 a32 a33 a34 b11 b12 b13 b14 b21 b22 b23 b24 b31 b32 b33 b34 A1=a11 a12 a13 a14 ;a21 a22 a23 a24; a31 a3
6、2 a33 a34, B1=b11 b12 b13 b14 ;b21 b22 b23 b24; b31 b32 b33 b34 运行,矩阵的运算(矩阵的加减、数乘、乘积),C=A1+B1 D=A1-B1 syms c, cA=c*A1 A2=A1(:,1:3), B1 G=A2*B1,矩阵的运算(矩阵的加减、数乘、乘积),A, A_trans=A H=1 2 3 ; 2 1 0 ; 1 2 3 , K=1 2 3 ; 2 1 0 ; 2 3 1 h_det=det(H), k_det=det(K), H_inv=inv(H),K_inv=K-1,矩阵的运算(左除和右除),左除“ ”：求矩阵方程
7、AX=B的解；（ A 、B的行要保持一致）解为 X=AB；当A为方阵且可逆时有X=AB=inv(A)*B； 右除“ / ”： 求矩阵方程XA=B的解 （A 、B的列要保持一致）解为 X=B/A ，当A为方阵且可逆时有X=B/A=B*inv(A),矩阵的运算(左除和右除),求矩阵方程： 设A、B满足关系式：AB2B+A,求B。 其中A=3 0 1; 1 1 0; 0 1 4。 解：有(A-2I)BA 程序 ： A=3 0 1; 1 1 0;0 1 4; B=inv(A-2*eye(3)*A, B=(A-2*eye(3)A 观察结果：,分块矩阵（矩阵的标识）,1.矩阵元素的标识 ：A(i,j)表示
8、矩阵A 的第 i 行 j 列的元素； 2.向量标识方式 A(vr,vc)：vr=i1,i2,ik、vc=j1,j2,ju分别是含有矩阵A的行号和列号的单调向量。A(vr,vc)是取出矩阵A的第i1,i2,ik行与j1,j2,ju列交叉处的元素所构成新矩阵。,分块矩阵（矩阵的标识）,取出A的1、3行和1、3列的交叉处元素构成新矩阵A1 程序 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分为四块，并把它们赋值到矩阵B中，观察运行后的结果。 程序 A11=A(1
9、: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，将它的第1行变为0。 程序： A=1 0 1 1 2;0 1 -1 2 3;3 0 5 1 0;2 3 1 2 1A(1,:)=0 0 0 0 0; A 删除上面矩阵A的第1、3行。 程序：A(1,3,:)= 结果,生成特殊矩阵,全1阵 ones(n), ones(m,n), ones(size(A) 全零阵：zeros(n) ，zeros(m,n), zeros(size(A) 常常用于对某个矩阵或向量
10、赋0初值 单位阵： eye(n)，eye(m,n) 随机阵： rand(m,n)， rand(n)=rand(n,n)用于随机模拟，常和rand(seed,k)配合使用。,生成特殊矩阵,将rand指令运行多次，观察结果。 程序： y1=rand(1,5), y2=rand(1,5), rand(seed,3), x1=rand(1,5), rand(seed,3), x2=rand(1,5) 结果,常用矩阵函数,det(A) ： 方阵的行列式； rank(A)： 矩阵的秩； eig(A)： 方阵的特征值和特征向量； trace(A)： 矩阵的迹； rref(A)： 初等变换阶梯化矩阵A svd
11、(A)： 矩阵奇异值分解。 cond(A)： 矩阵的条件数；,数据的简单分析,1.当数据为行向量或列向量时，函数对整个向量进行计算. 2.当数据为矩阵时，命令对列进行计算，即把每一列数据当成同一变量的不同观察值。 max(求最大)、min(求最小)、mean(求平均值)、sum(求和)、std(求标准差)、cumsum(求累积和)、median(求中值)、diff(差分)、sort(升序排列)、sortrows(行升序排列)等等。,数据的简单分析,观察：生成一个36的随机数矩阵，并将其各列排序、求各列的最大值与各列元素之和。 程序 rand(seed,1);A=rand(3,6), Asort
12、=sort(A), Amax=max(A), Asum=sum(A) 结果,实验二,函数可视化与Matlab作图,函数的可视化,f (x), g (x)是周期函数吗？观察它们的图象。程序 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),令,绘制平面曲线(plot指令),plot
13、(x,y)： 以x为横坐标、y为纵坐标绘制二维图形 x,y是同维数的向量； plot(y)： 相当于x=1,2,length(y)时情形。,绘制平面曲线（绘制多个图形),1. plot(x,y1;y2;),x是横坐标向量，y1;y2;是由若干函数的纵坐标拼成的矩阵 2. plot(x,y1), hold on, plot(x,y2), hold off 3. plot(x,y1,x,y2,) 4.plotyy 两个坐标系，用于绘制不同尺度的函数。,绘制平面曲线（线型、点形和颜色的控制）,plot(x,y,颜色线型点形) plot(x,y,颜色线型点形,x,y,颜色线型点形, ) 句柄图形和se
14、t命令改变属性值，可套用： h=plot(x,y),set(h,属性,属性值,属性,属性值,) 也可用plot(x,y,属性,属性值)设置图形对象的属性。,绘制平面曲线（属性变量和属性值）,线宽：LineWidth 点的大小： MarkerSize 线型：LineStyle 颜色：color,绘制平面曲线（例）,观察： 改变绘图的线型和颜色。 用grid on 指令为图形窗口加上 网格线，并改变网格的线型和字体的大小。 程序 h=plot(0:0.1:2*pi,sin(0:0.1:2*pi); set(h,LineWidth,5,color,red); grid onset(gca,GridL
15、ineStyle,-,fontsize,16) 观察结果,绘制平面曲线（坐标轴的控制）,axis指令axis(xmin xmax ymin ymax)：设定二维图形的x和y坐标的范围；axis(xmin xmax ymin ymax zmin ymax)：设定三维图形的坐标范围 ； 其中xminxxmax， yminyymax ，zminzzmax。,绘制平面曲线（gca属性控制）,改变当前轴对象句柄gca属性 用set(gca,属性,属性值,)可改变字体大小、坐标刻度等轴对象的内容。例如： set(gca,ytick,-1 -0.5 0 0.5 1) 将 y 坐标按向量-1 -0.5 0 0
16、.5 1将刻度分成4格； set(gca,yticklabel,a|b|c|d|e) 改变y坐标刻度的说明。,绘制平面曲线（gca属性控制,例）,设置y坐标的刻度并加以说明，并改变字体的大小。 程序 plot(0:0.1:2*pi,sin(0:0.1:2*pi),k.-,);grid on,axis(0 6.3 -1.1 1.1)，set(gca,ytick,-1 -0.5 0 0.5 1), set(gca,yticklabel,a|b|c|d|e),set(gca,fontsize,20)get(gca) 运行结果,绘制平面曲线（文字标注）,title(图形标题)； xlabel（x轴名称
17、）；ylabel（y轴名称）；zlabel（z轴名称）； text(说明文字)：创建说明文字； gtext（说明文字）：用鼠标在特定位置输入文字。 文字标注常用符号：pi （）；alpha （）；beta （）； leftarrow （左箭头） rightarrow （右箭头）；bullet （点号）,绘制平面曲线（程序讲解，exp2_1.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,ti
18、mes New Roman), xlabel(itt(deg); ylabel(itmagnitude); title( itsine wave and itAe-alphaittwave);,绘制平面曲线（程序讲解，exp2_1.m）,text(6,sin(6),fontsize15The Value itsin(t) at itt=6rightarrowbullet, HorizontalAlignment,right), text(2,3*exp(-0.5*2), fontsize15bulletleftarrow The Value of it3e-0.5 itt=, num2str(
19、3*exp(-0.5*2), at itt =2 ); legend(itsin(t),itAe-alphat) 注1： num2str： string1 ,num2str,string2，用方括号 注2： legend 请结合图形观察此命令的使用,图形窗口的创建和分割,subplot(m,n,k)命令。在图形区域中显示多个图形窗口。m为上下分割数，n为左右分割数，k为第k子图编号。 例：将一个图形分为9个子图，在第k个子图画sin(kx) 的图象. 程序： clf，b=2*pi;x=linspace(0,b,50);for k =1:9y=sin(k * x);subplot(3,3,k),
20、plot(x,y),axis(0,2*pi,-1,1) end,若干有用的指令,clf：清除图形窗口已有的内容. shg：显示图形窗口。 clear、 clear x：清除工作空间的已有变量。 figure(n): 打开第n个图形窗口 help: : 续行号,绘制二元函数,基本步骤： 1.生成二维网格点 2. 计算函数在网格点上的值 3. 绘制函数图形,三维绘图（ meshgrid指令：生成网格点）,观察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=meshg
21、rid(x,y); plot(X,Y,+) 观察结果,三维绘图（计算函数值，定义域裁减）,程序： for i=1:n for j=1:nif (1-X(i,j)eps1|X(i,j)-Y(i,j)eps1 z(i,j)=NaN; elsez(i,j)=1000*sqrt(1-X(i,j)-1.*log(X(i,j)-Y(i,j); endend end,三维绘图（绘图指令）,mesh(X,Y,z) ： 在三维空间中绘出由(X,Y,z)表示的曲面； meshz(X,Y,z)： 除了具有mesh的功能外，还画出上下高度线， meshc(X,Y,z)： 除了具有mesh的功能外，还在曲面的下方画出函
22、数z=f(x,y)的等值线图， surf(X,Y,z)： 也是三维绘图指令，与mesh的区别在于mesh绘出彩色的线，surf绘出彩色的面， 运行exp2_1,观察效果,三维绘图（等值线指令）,表现二维函数的图形的另一种方式是绘制等值线图。 contour(X,Y,z,n)： n条等高线，n可缺省； contourf(X,Y,z,n)： 等值线间用不同的颜色填满，有更好的视觉效果； contour3(X,Y,z,n)： 在三维空间画出等值线图 colorbar： 将颜色与函数值对应起来显示在图中。,三维绘图（等值线指令，继续exp2_2显示效果）,clf,contour(X,Y,z,40),c
23、olorbar contourf(X,Y,z,40),colorbar contour3(X,Y,z,40),colormap(0,0,0) 为等值线标上函数值: 可套用下面程序的格式.cs,h=contour(X,Y,z,15); clabel(cs,h,labelspacing,244) labelspace是数值标记之间相隔的宽度，默认值为144, 这里取了244，,空间曲线和运动方向的表现,一条空间曲线可以用矢量函数表示为,它的速度矢量表现为曲线的切矢量:,空间曲线和运动方向的表现,很显然飞行曲线方程为：,绘制空间曲线（指令）,plot3(x,y,z)： 绘制三维空间曲线，用法和plo
24、t类似。 quiver(X,Y,u,v)：绘制二维矢量， 在坐标矩阵点X,Y处绘制矢量u,v, 其中u为矢量的x坐标，v为矢量的y 坐标，其维数不小于2。 quiver3（X,Y,Z,u,v,w）： 绘制三维矢量，用法与quiver类似。 Gradient：Fx,Fy,Fz=gradient(F)为函数F数值梯度,绘制空间曲线（程序讲解exp2_3）,exp2_3.m clf,t=linspace(0,1.5,20); x=t.2;y=(2/3)*t.3;z=(6/4)*t.4-(1/3)*t.3; plot3(x,y,z,r.-,linewidth,1,markersize,10),hold
25、 on Vx=gradient(x);Vy=gradient(y);Vz=gradient(z); h=quiver3(x,y,z,Vx,Vy,Vz),set(h,linewidth,1),grid on axis(0 1.5 0 1.5 0 40) xlabel(x),ylabel(y),zlabel(z)，box on 运行程序,应用、思考和练习,用平行截面法讨论由曲面z=x2-y2构成的马鞍面形状。 对于二重积分，积分指示线方法是很有用的，当然你首先得了解一下什么是积分指示线法，请查阅高等数学相关的内容，然后设计一个数学实验，然后用Matlab的绘图工具表现这一方法。,应用、思考和练习,
26、绘制微分方程 y/dx=xy, y(0)=0.4的斜率场， 并将解曲线画在图中，观察斜率场和解曲线的关系。,应用、思考和练习,地球表面的气温差异很大，而且随时间变化，要绘制气温分布图绝不是一件容易的事情。但是赤道温度最高，而在两级最冷，在中间地带则是过渡带。所以可粗略将这种气温分布情况用图形表现出来，试绘制地球表面的气温分布示意图。,应用、思考和练习,应用、思考和练习(若干特殊图形 ),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 1
27、0 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),应用、思考和练习(周期函数的推广）,应用、思考和练习(线性P周期函数）,应用、思考和练习(线性P周期函数),应用、思考和练习(线性P周期函数）,p周期函
28、数是一种特殊的线性p周期函数，为M=0的情形。 线性p周期函数是不是一个线性函数同一个 p 周期函数之和呢？ 如何证明这一点？,应用、思考和练习（循环比赛的名次）,有若干支球队参加循环比赛，他们两两交锋，每场比赛只计胜负，不允许平局，循环赛结束后要根据他们的比赛成绩排列名次。 一种方法是计算得分，得分是每支球队获胜的场次，根据各队的得分排出名次，决定冠军队。,应用、思考和练习（循环比赛的名次,虽然计算各队的得分很容易，但有时按得分排名的方法并不一定合理。 假定有4支球队, 记做v1v4。在一次循环赛中，v1得分为2，v2得分为2，v3得分为1，v4得分为1，可以把得分写成4维向量的形式：s =
29、 2 2 1 1T， 在这种情况下应该如何决定v1, v2,v3,v4的名次呢？,函数式直接确定型模型,实验三,从系统分析的观点理解函数,y = f（x）x：自变量，y：因变量，f： 映射规则,函数不是枯燥的数学符号 函数是系统 函数是数学模型 是描述自然现象的有力工具,黑箱模型和经验函数,白箱：映射规则f 已知； 灰箱：映射规则f 部分已知； 黑箱：映射规则f 未知。 对于黑箱模型，只知道输入输出的数据， 需根据这些数据近似决定映射规则f,经验函数(机床加工问题）,用程控铣床加工机翼断面的下轮廓线时 每一刀只能沿x方向和y方向走非常小的一步。 表3-1给出了下轮廓线上的部分数据 但工艺要求铣
30、床沿x方向每次只能移动0.1单位. 这时需求出当x坐标每改变0.1单位时的y坐标。 试完成加工所需的数据，画出曲线.,交通事故的调查（司机有责任吗？）,一辆汽车在拐弯时急刹车，结果冲进路边的沟里 警察闻讯赶到现场，对汽车留在路上的刹车痕迹进行了测量，确定刹车痕迹近似为一圆弧。 讯问司机时，他说，当车进入弯道时刹车失灵，且进入弯道时汽车时速为40英里小时。 此速度系该路的速度上限。 通过验车证实该车制动器在事故发生时确实失灵， 但司机所说的车速是否真实呢？,交通事故的调查（司机有责任吗？）,通常，作一条基准线来测量刹车痕迹. （水平）距离x沿基准线测得，（垂直）距离y与x垂直. 表3-6给出了刹
31、车痕迹的测量有关值. 警察测量了路的坡度，发现这段路是平的 请给出一个使警察可以核对速度的计算办法,航行区域的警示线,某海域上频繁地有各种吨位的船只经过。 为保证船只的航行安全，有关机构在低潮时对水深进行了测量,表3-8是他们提供的测量数据： 表3-8. 水道水深的测量数据 x 129.0 140.0 103.5 88.0 185.5 195.0 105.5 y 7.5 141.5 23.0 147.0 22.5 137.5 85.5 z 4 8 6 8 6 8 8 x 157.5 107.5 77.0 81.0 162.0 162.0 117.5 y -6.5 -81.0 3.0 56.5
32、-66.5 84.0 -33.5 z 9 9 8 8 9 4 9,航行区域的警示线,其中（x, y）为测量点，z为（x, y）处的水深(英尺)。 船的吨位可以用其吃水深度来反映， 分为 4英尺、4.5英尺、5英尺和 5.5英尺 4 档。航运部门要在矩形海域（75，200）（50，150）上为不同吨位的航船设置警示标记。 请根据测量的数据描述该海域的地貌，并绘制不同吨位的警示线，供航运部门使用。 水深z是区域坐标（x, y）的函数z= z (x, y）， 测量数据只是它的部分取值。 可绘制函数图象和等值线图，将不同吃水线标记图上,插值与拟合（基本原理和区别）,已知有n +1个节点（xj,yj），
33、j = 0,1, n 其中xj互不相同 节点(xj, yj)可看成由某个函数 y= f（x）产生 f 的解析表达式可能十分复杂 或不存在封闭形式, 也可以是未知的,插值与拟合（基本原理和区别）,插值： 构造一个相对简单的函数 y=g(x) 使g通过全部节点 即使g (xj) = yj，j=0,1, n 用g (x)作为函数f ( x )的近似。 插值指令 yi=interp1(x1,y1,xi,method) 对应于插值函数yi=g(xi) 其中x1,y1为节点向量method四个选项： nearest 为近邻插值；linear为线性插值； spline 为样条插值； cubic为立方函数插值
34、。,插值与拟合（基本原理和区别）,多项式拟合 对给定的数据（xj,yj），j = 0,1, n 选取适当阶数的多项式(也可采用其它形式的函数) 例如二次多项式g(x)=ax2+bx+c 使g(x)尽可能逼近（拟合）这些数据 拟合指令polyfit、polyval 用p=polyfit(x1,y1,m)做 m 次多项式拟合 拟合数据向量为x1,y1 多项式系数为p=p(1),p (m),p (m+1) 即g(x)=p(1)xm+p (m)x+p (m+1) 用y = polyval(p,x)计算在x 处 多项式的值 y,观察插值、拟合的效果,运行观察程序exp3_1.m 选取一个已知函数作为参考
35、，并将这一函数的图象用虚线显示在图中。 观察程序允许用鼠标选取节点 按鼠标左键选点，按右键选最后一个点 观察不同的选点方式对各种插值和拟合效果的影响,程序注解（inline指令）,定义内联函数：inline指令 g=inline(x2-x4); 程序,程序注解（ginput）,x,y,button = ginput(n) 用鼠标在屏幕选n个点，返回这n个点，存于x,y中。 button 记录了选点时使用的鼠标键方式： 1为左键、2为中间键、 3为右键。,程序注解（插值拟合）,xx=linspace(a,b,n); %定义自变量xx ynearest=interp1(x1,y1,xx,neare
36、st);ylinear=interp1(x1,y1,xx,linear); yspline=interp1(x1,y1,xx,spline);p,c=polyfit(x1,y1,4); ypolyfit=polyval(p,xx);,程序注解（插值拟合）,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)
37、 subplot(2,2,4), h=plot(xx,ypolyfit,r-);set(h,linewidth,2),插值拟合效果观察,沿曲线选取个节点，保持等间隔。当节点较少时，插值的效果如何？ 加密节点，共个等距节点，观察插值的效果，如果去掉中间的一个节点，插值效果又会如何？ 有意偏离原来的曲线，如果误差较大，将会怎样呢？,微分、积分和微分方程,实验四,定积分-连续求和,定积分-连续求和,三种方法计算数值积分,（1）定义法，取近似和的极限。 高等数学中不是重点内容 但数值积分的各种算法却是基于定义建立的（2）用不定积分计算定积分。 不定积分是求导的逆运算， 而定积分是连续变量的求和（曲边梯
38、形的面积） 表面上看是两个完全不同的概念， 通过牛顿莱布尼兹公式联系在一起， （3）解微分方程计算定积分,微积分学基本定理,特别，F(b)-F(a) 就是所需的定积分. 在高等数学中总是期望求出不定积分的封闭解. 但数值积分是更有用的工具。 牛顿莱布尼兹公式不愧为微积分的“基本定理”。,基本定理的推广（解微分方程计算定积分）,基本定理的推广（解微分方程计算定积分）,解微分方程的 Eular折线法,解微分方程的 Eular折线法,将区间n = 4等分（共有5个分点)；计算分点和相应的函数值（x(1),x(2), x(3) x(4) x(5)（f(1),f(2) ,f(3) ,f(4) ,f(5)
39、 在第一个子区间x(1),x(2)上，画出折线段 y(2)=y(1)+f(1)*(x-x(1)代替解曲线段y(x)， 这里y(1)=y0=0 折线段的起点为x(1), y(1)，终点为x(2),y(2). 运行exp4_1.m，观察第二、三、四子区间的情况。,符号微积分,用Matlab符号工具箱（Symbolic Toolbox）可以进行符号演算,符号微积分(创建符号变量）,sym var 创建单个符号变量；syms var1 var2 创建多个符号变量； f=sym(符号表达式) 创建符号表达式，赋予f； equ=sym(equation) 创建符号方程 。,符号微积分(极限）,limit(
40、表达式,var,a)：求当var a，表达式的极限 例：求极限：,syms x a I1=limit(sin(x)-sin(3*x)/sin(x),x,0) 运行结果,符号微积分(求导）,diff(f,var,n) 求 f 对变量var 的n阶导数 缺省n时为求一阶导数 缺省变量var 时，默认变量为x 可用来求单变量函数导数 多变量函数的偏导数 还可以求抽象函数的导数,符号微积分(求导）,例：求,syms x y f=sym(exp(-2*x) * cos(3 * x(1/2) diff(f,x) 运行,符号微积分(求导）,syms x y g=sym(g(x,y) f=sym(f(x,y,
41、g(x,y) diff(f,x) diff(f,x,2) 运行,例：求,符号微积分(积分）,int(f,var)：求函数f的不定积分； int(f, var, 积分下限，积分上限)： 求函数f的定积分或广义积分 例：求不定积分,syms x y zI1=int(sin(x*y+z),z),符号微积分(积分）,syms x y zI2=int(1/(3+2*x+x2),x,0,1)I3=int(1/(3+2*x+x2),x,-inf,inf),符号微积分(化简、提取和代入),符号运算的结果比较繁琐，使用化简指令可对其进行化简。 但是不能指望机器可以完成一切，人的推理往往必须的。 常用的化简指令如
42、下展开指令：expand(表达式)；因式分解：factor(表达式)降幂排列：collect（表达式，var） ； 一般化简：simplify(A)；,符号微积分(化简、提取和代入),观察： 将展开(a+x)6-(a-x)6，然后作因式分解。 t_expand=expand(t) t_factor=factor(t_expand) t_simplify=simplify(t) 观察结果,数值微积分（梯形公式和辛普森公式）,trapz(x,y)，按梯形公式计算近似积分； 其中步长x=x0 x1 xn和函数值y=f0 f1 fn为同维向量， q = quad(fun,a,b,tol,trace,P
43、1,P2,.) （低阶方法，辛普森自适应递归法求积） q = quad8(fun,a,b,tol,trace,P1,P2,.)（高阶方法，自适应法Cotes求积） 在同样的精度下高阶方法quad8要求的节点较少。,x,y=ode23(fun,tspan,y0,option) （ 低阶龙格库塔函数） x,y=ode45(fun,tspan,y0,option) （高阶龙格库塔函数）,应用、思考和练习(追击问题),我缉私雷达发现，距离d处有一走私船正以匀速a沿直线行驶，缉私舰立即以最大速度（匀速v）追赶。 若用雷达进行跟踪，保持船的瞬时速度方向始终指向走私船， 缉私舰的运动轨迹是怎样的？是否能够追
44、上走私船？ 如果能追上，需要用多长时间？,应用、思考和练习(追击问题),应用、思考和练习(追击问题),r = dsolve(eq1,eq2, cond1,cond2, v) 方程的符号解 syms y d r xs1= dsolve(D2x=-r*sqrt(1+Dx2)/y,x(20)=0,Dx(20)=0,y)xs=simplify(xs1) 运行结果，画彗星图,应用、思考和练习(追击问题),r = dsolve(eq1,eq2, cond1,cond2, v) 方程的符号解 syms y d r xs1= dsolve(D2x=-r*sqrt(1+Dx2)/y,x(20)=0,Dx(20)
45、=0,y)xs=simplify(xs1) 运行结果，画彗星图,应用、思考和练习(追击问题，如果雷达失效),当缉私舰雷达发现d处有一走私船后，雷达突然损坏 若假定走私船作匀速直线运动（但不知方向），且缉私舰艇速度v大于走私船速度a， 则缉私舰应采用什么样的航行路线，不管走私船从哪个方向逃跑，都能追捕上它？,实时动画制作（见实验10）,观察：模拟弹簧振动 讨论最简单的情形，一弹簧系统作横向运动，其位移由u=2+cos(t) 所决定， 仿真弹簧的振动,实时动画制作(初始化、见实验10）,程序讲解 animinit(onecart1 Animation) axis(-2 6 -10 10); hol
46、d on; u=2;xy= 0 0 0 0 u u u+1 u+1 u u;-1.2 0 1.2 0 0 1.2 1.2 -1.2 -1.2 0;x=xy(1,:);y=xy(2,:); plot(-10 20,-1.4 -1.4,k-,LineWidth,2); hndl=plot(x,y,k-,EraseMode,XOR,LineWidth,2),实时动画制作(初始化、见实验10）zxy10-2,set(gca,UserData,hndl); for t=1:0.025:1000;u=2+exp(-0.00*t)*cos(t);x=0 0 0 0 u u u+1 u+1 u u;hndl=
47、get(gca,UserData);set(hndl,XData,x,YData,y);drawnow end,电影动画制作（zxy7_3),moviein、 getframe、movie指令 x=-8:0.5:8; XX,YY=meshgrid(x); r=sqrt(XX.2+YY.2)+eps; Z=sin(r)./r; surf(Z); %画出祯 theAxes=axis; %保存坐标值，使得所有帧都在同一坐标系中,电影动画制作,fmat=moviein(20); %创建动画矩阵，保存20祯 for j=1:20; %循环创建动画数据 surf(sin(2*pi*j/20)*Z,Z) %画出每一步的曲面axis(theAxes) %使用相同的坐标系 fmat(:,j)=getframe;%拷贝祯到矩阵fmat中 end movie(fmat,10) %演示动画10次,
