第3章_Matlab程序设计
第三章Matlab 绘图与例题 (1)
[X ,Y ] meshgrid(x, y)
mesh(X ,Y , Z )
• mesh(Z) 若提供参数x,y,等价于mesh(x,y,Z),否 则默认
• mesxhc 1除: n了, y生1成: m网格曲面，还在xoy面上生成等高线图
• meshz 除了生成网格曲面,还在曲线下加上一个长方形 台柱。
机动 目录 上页 下页 返回 结束
空间网格曲面。而函数meshgrid就用来生成xoy 平面上的小矩形顶点坐标值的矩阵（格点矩阵）。 函数meshgrid也适用于三元函数。
meshgrid的调用形式
• [X,Y]=meshgrid(x,y) 绘制二维图形时生成小矩形的格点 • [X,Y]= meshgrid(x) 等价于[X,Y]=meshgrid(x,x) • [X,Y,Z]=meshgrid(x,y,z)绘制三维图形生成小矩形的格点 • [X,Y,Z]= meshgrid(x) 等价于[X,Y,Z]=meshgrid(x,x,x)
c青
x叉
> 右三角形 -- 划线
m 洋红 + 十字 < 左三角形
y黄
*星
p 五角形
k黑
s 方块 h 六角形
d 菱形
例、在同一图形窗口画出不同颜色，线形的正、余 弦图形
机动 目录 上页 下页 返回 结束
三、 基本绘图控制参数
常用作图指令 • 清除图形窗口 clf • 显示/不显示格栅 grid on/off • 保留/释放现有图形 hold on/off • 使用/不使用坐标框 box on/off • 暂停直到击键盘 pause
例3-1-1、在同一坐标系中画出两个函数 y cos 2x, y x 的图形，自变量范围为：2 x 2 ,及5个点
matlab第一讲
·
1.2 MATLAB 集成开发环境
如果一个命令行很长，一个物理行之内写不下，可以在第1个物理行之 后加上3个小黑点并按下回车键，然后接着下一个物理行继续写命令的 其他部分。3个小黑点称为续行符，即把下面的物理行看作该行的“逻 辑”继续。 例如：
z=1+1/(1*2)+(1*2*3)+1/(1*2*3*4)+…
系统科学研究所
1.2 MATLAB 集成开发环境
MATLAB的工具栏提供了一些命令按钮和一个当前路径列表框。
2.命令窗口
一般来说，在命令编辑区的一个命令行输入一条命令，命令行以回车 结束。但一个命令行可以输入若干条命令，各命令之间以逗号分隔， 若前一命令后带有分号，则逗号可以省略。例如：
x=720,y=68 x= 720 y=86 x=720;y=86 y=86
1.3MATLAB 的帮助功能
1.3.1帮助界面
进入MATLAB帮助界面可以通过以下3种方法。 单击MATLAB主窗口工具栏中的Help按钮。
在命令窗口中输入”helpwin”、”helpdesk”或“doc”命令。
选择Help菜单中的”MATLAB Help”选项。 1.3.2帮助命令
1.help命令
系统科学研究所
1.2 MATLAB 集成开发环境
检查该命令是否为MATLAB 搜索路径中其他目录下的M文件。 2.设置搜索路径
用户可以将自己的工作目录列入MATLAB搜索路径，从而将用户目录
纳入MATLAB系统统一管理。 （1）用path命令设置搜索路径 （2）用对话框设置搜索路径
·
系统科学研究所
统程序设计语言一样进行程序设计，而且结合MATLAB的数值计算和
第三章_matlab矩阵运算
主讲：陈孝敬 E-mail:chenxj9@
第3章
数学运算
主要内容：
①矩阵运算； ②矩阵元素运算；
3.1 矩阵运算
3.1.1 矩阵分析
1．向量范式定义：
x x x
1
n
k 1
xk
2 k
2
k 1 n
x
n
1/ 2
k 1
xk
向量的3种常用范数及其计算函数 在MATLAB中，求向量范数的函数为： (1) norm(V)或norm(V,2)：计算向量V的2—范数。 (2) norm(V,1)：计算向量V的1—范数。 (3) norm(V,inf)：计算向量V的∞—范数。
3.1.2 矩阵分解
矩阵分解：把矩阵分解成比较简单或对它性质比较熟悉的若干 矩阵的乘积的形式；
1．Cholesky分解： Cholesky分解是把对称正定矩阵表示成上三角矩阵的转 置与其本身的乘积，即：A=RTR，在Matlab中用函数chol 来计算Cholesky分解 例3-13 求矩阵A=pascal（4）的Cholesky分解， A=pascal(4) R=chol(A) R’*R
例3-18.求解方程组
x1 x2 3 x3 x4 1 3 x1 x2 3 x3 4 x4 4 x 5x 9 x 8x 0 2 3 4 1
解 先用Matlab函数null求出对应的齐次线性方程组的基础解 系,再利用其系数矩阵的上、下三角阵求出方程组的一个特解, 这样即可得到该方程组的通解,程序如下: >> >> >> >> >> >> A=[1 1 -3 -1;3 -1 -3 4;1 5 -9 -8]; b=[1 4 0] ′; format rat C=null(A , ′r′); %求基础解系 [L,U]=lu(A); %A=LU,L为上三角阵,U为下三角阵 X0= U\(L\b) %用LU求出一个齐次方程的特解
matlab电子教案
第1章MATLAB操作基础1.1 MATLAB概述1.2 MATLAB的运行环境与安装1.3 MATLAB集成环境1.4 MATLAB帮助系统1.1 MATLAB概述•1.1.1 MATLAB的发展1984年，MATLAB第1版(DOS版)1992年，MATLAB4.0版1994年，MATLAB 4.2版1997年，MATLAB 5.0版1999年，MATLAB 5.3版2000年，MATLAB 6.0版2001年，MATLAB 6.1版2002年，MATLAB 6.5版2004年，MATLAB 7.0版•1.1.2 MATLAB的主要功能1．数值计算和符号计算功能MATLAB以矩阵作为数据操作的基本单位，还提供了十分丰富的数值计算函数。
MATLAB和著名的符号计算语言Maple相结合，使得MATLAB具有符号计算功能。
2．绘图功能MATLAB提供了两个层次的绘图操作：一种是对图形句柄进行的低层绘图操作，另一种是建立在低层绘图操作之上的高层绘图操作。
•3．编程语言MATLAB具有程序结构控制、函数调用、数据结构、输入输出、面向对象等程序语言特征，而且简单易学、编程效率高。
4．MATLAB工具箱MATLAB包含两部分内容：基本部分和各种可选的工具箱。
MATLAB工具箱分为两大类：功能性工具箱和学科性工具箱。
•1.1.3 初识MATLAB例1-1 绘制正弦曲线和余弦曲线。
x=[0:0.5:360]*pi/180;plot(x,sin(x),x,cos(x));•例1-2 求方程3x4+7x3 +9x2-23=0的全部根。
p=[3,7,9,0,-23]; %建立多项式系数向量x=roots(p) %求根•例1-3 求积分quad('x.*log(1+x)',0,1)•例1-4 求解线性方程组。
a=[2,-3,1;8,3,2;45,1,-9];b=[4;2;17];x=inv(a)*b1.2 MATLAB的运行环境与安装•1.2.1 MATLAB的运行环境硬件环境：(1) CPU(2) 内存(3) 硬盘(4) CD-ROM驱动器和鼠标。
常用数学软件教程 033 第3章 MATLAB使用基础 第3节 字符串和符号矩阵
3.3字符串和符号矩阵在数学、物理等许多基础学科和电路、力学工程应用中，常常要进行公式的推演变换、代数方程的符号解。
微积分的解析解等，处理这类问题就要用到符号运算。
符号运算是对方程或代数式中的数学符号进行运算操作，运算结果是数学表达式而不是数值。
例如高等数学中有求函数的导函数和求函数导数在某一点取值的问题，前者属于符号运算，后者属于数值计算。
MA TLAB开发了符号运算工具包（Symbolic Math Toolbox），可以进行代数式的因式分解。
展开和化简，函数的幂级数展开，求微积分及微分方程组的解析解等多种运算。
另外，MA TLAB中的符号运算是按有理数的运算方式进行的，如对未知量x、有理数（如1/3）或超越数（如π）等运算都没有有效位数的限制，也没有舍入误差，所以相当精确，这种精确的结果可以任意转换成任意精度的数值解，这对数值计算来说，是一个难能可贵的优点。
3.3.1字符串和自定义函数字符串在数据处理、构造表和函数求值运算中都很有用。
为了按照要求的格式进行结果的显示，就需要熟悉MA TLAB中关于字符串的操作。
1.字符串及其显示（1）字符串和字符串变量名。
用单引号界定的排成序列的各种符号，包括数字、英文、横线、括号及表达式、方程式等，都被看作字符串（也称字符串数据或支付量）。
用赋值号“=”据按一定的格式显示。
a）字符串数据显示指令。
最常用的指令是disp，它的使用格式为：●disp(zs)①输入参数zs可以是数值变量、字符串变量（或它们的变量名）和被界定的字符串，但每次只能有一个内容，各变量名之间用空格分开。
②zs若是字符串，回车后照原样显示（包括字符间空格大小）；zs若是数值变量名，则显示它代表的数据。
b）空字符（空格）输出显示指令。
字符或数据间的空格大小，除在字符串间留有空格外，也可以利用指令blanks控制，它的使用格式为：●blanks (n)①该指令作为disp指令的输入参数，可夹杂在数据或字符串中间使用。
MATLAB编程
第一章入门MATLAB的灵活性和平台独立性是通过将MATLAB代码编译成设备独立的P代码，然后在运行时解释P代码来实现的。
但MATLAB是解释型语言不是编译型语言，产生的程序执行速度慢。
MATLAB有两个基本的缺点：1.是解释型语言，其执行速度要比编译型语言慢得多。
这个问题可以通过合理的MATLAB结构得到缓解，也可以在发行广泛使用前编译出MATLAB程序2.费用高MATLAB的开发环境任何一个MATLAB程序的基本组成单元是数组。
数组是一组数据值得集合，这些数据被遍上行号和列号，拥有唯一的名称。
数组中的单个数据是可以通过带有小括号的数组名访问的，括号内有这个数据的行标和列标，中间用逗号隔开。
标量也被MATLAB当作数组，只不过只有一行和一列。
如果一个语句在一行内书写太长了，可能要另起一行接着写，这种情况需要在第一行末打上半个省略号（…），再开始第二行的书写。
历史命令窗口（The history command window）顺序是按逆序排列的，双击这些命令可以再次执行。
启动平台（The launch pad）是一个特殊工具，为MATLAB和其工具箱提供帮助、demo、其他相关文件和应用程序等参考资料。
编辑调试器，一般用于创造M文件或者修改已存在的M文件，当打开或修改一个M文件编辑调试器会被自动调用。
M文件中的评论用绿色表示，变量和数字用黑色来表示，字符变量用紫色表示，语言的关键字用蓝色表示。
当一个专门的命令、M文件或函数运行时，工作区是MATLAB所需要的所有变量和数组的集合。
所有在命令窗口中执行的命令和所有在命令窗口执行的脚本文件（Script files）都会被分配一个普通的分配空间，所以它们能共享变量。
MATLAB 函数拥有独立的工作区，这是函数区别于脚本文件的一个重要特征。
Who命令会产生一个在当前工作区内的所有变量和数组状况表一些重要命令Clc、clf、clear！的特殊作用是给计算机操作系统发送一个命令。
杭州电子科大,MATLAB经典教程。第三章 数据类型
例子3-13 复杂的关系运算。 在MATLAB命令行窗口中，键入下面的指令： >> A = reshape(-4:4,3,3) A= -4 -3 -1 0 2 3
-2
1
4
>> B = ~(A>=0) （提问：输出的结果是什么？）
3.3.4 运算符的优先级
M语言的运算符和普通的高级编程语言类似，也具有 相应的计算优先级：
note： 1) eps、realmax和realmin三个常量具体的数值与运行 MATLAB的计算机相关，不同的计算机系统可能具有不同 的数值。 2) MATLAB的常量数值是可以修改的。例如，pi=100， 这样pi这个常量的数值就变成了100。但是，如果用clear指 令清除变量后，pi将变成原有的常量数值。
例子3-10 逻辑运算示例。
在MATLAB命令行窗口中，键入下面的指令：
>> a = eye(3);
>> b = a;b(3,1) =1;
>> a && b (标量)
??? Operands to the || and && operators must be
convertible to logical scalar values. >> a & b 例子3-11 函数all和any使用示例。
例子3-2 使用不同的数据类型。 在MATLAB命令行窗口中，键入下面的指令： >> a = [ 1 2 3]; >> b = [ 3 4 5]; >> c = a+b; >> int16(a)+int16(b) ??? Error using ==> + Function '+' is not defined for values of class 'int16'. note：MATLAB系统默认的运算都是针对双精度类型，所 以在进行两个int16类型的变量加法时，MATLAB报告了相 应的错误。（解决方法：计算前后进行类型转化。）
汽车理论课后习题答案及MATLAB编程
汽车理论课后习题答案—附MATLAB编程教材：汽车理论（第5版）清华大学余志生主编机械工业出版社出版本文档包含第三章 4.3 5.11 6.2 6.4 6.5所要求的MATLAB编程其它习题答案：第一章：1.1 1.2 1.7 1.8第二章：2.2 2.3 2.4第三章第四章：4.1 4.2 4.3 4.5第五章：5.1 5.8 5.10 5.11 5.12 5.14 5.17第六章：6.1~6.5见本作者另一百度文库文档。
获取方法：点击左上角本作者的昵称，进入本作者的文库，《汽车理论课后习题答案（含MATLAB编程）》第三章：m=3880;%总质量g=9.8;%重力加速度r=0.367;%车轮半径eta_t=0.85;%传动系机械效率f=0.013;%滚动阻力系数CDA=2.77;%空气阻力系数*迎风面积i0=[5.175.435.836.176.33];%主减速器传动比If=0.218;%飞轮转动惯量Iw1=1.798;Iw2=3.598;%两前轮/四后轮转动惯量Iw=Iw1+Iw2;ig4=[6.09 3.091.711];%4挡变速器传动比ig5=[5.56 2.7691.64410.793];%5挡变速器传动比n=600:1:4000;%发动机转速Tq=-19.313+295.27*n/1000-165.44*(n/1000).^2+40.874*(n/1000).^3-3.8445*(n/1000).^4;%4挡发动机for i=1:5Ft1=Tq*ig4(1)*i0(i)/r;%各档对应转速下的驱动力Ft2=Tq*ig4(2)*i0(i)/r;Ft3=Tq*ig4(3)*i0(i)/r;Ft4=Tq*ig4(4)*i0(i)/r;u1=0.377*r*n/(i0(i)*ig4(1));%各档对应转速下的车速u2=0.377*r*n/(i0(i)*ig4(2));u3=0.377*r*n/(i0(i)*ig4(3));u4=0.377*r*n/(i0(i)*ig4(4));F1=m*g*f+CDA*u1.^2/21.15;%各档对应转速下的行驶阻力F2=m*g*f+CDA*u2.^2/21.15;F3=m*g*f+CDA*u3.^2/21.15;F4=m*g*f+CDA*u4.^2/21.15;delta1=1+(Iw+If.*ig4(1).^2.*i0(i).^2.*eta_t)./(m.*r.^2);%汽车旋转质量换算系数delta2=1+(Iw+If.*ig4(2).^2.*i0(i).^2.*eta_t)./(m.*r.^2);delta3=1+(Iw+If.*ig4(3).^2.*i0(i).^2.*eta_t)./(m.*r.^2);delta4=1+(Iw+If.*ig4(4).^2.*i0(i).^2.*eta_t)./(m.*r.^2);a1=(Ft1-F1)/(delta1*m);%加速度a2=(Ft2-F2)/(delta2*m);a3=(Ft3-F3)/(delta3*m);a4=(Ft4-F4)/(delta4*m);b1=1./a1;%加速度倒数b2=1./a2;b3=1./a3;b4=1./a4;figure(1)subplot(2,3,i)plot(u1,b1,u2,b2,u3,b3,u4,b4);title('加速度倒数-速度曲线');xlabel('u')ylabel('1/a')gtext('1/a1')gtext('1/a2')gtext('1/a3')gtext('1/a4')gtext(['i0='num2str(i0(i))]);%计算采用2挡起步加速至96.6km/h的原地起步加速时间u1min=min(u1);u1max=max(u1);u2min=u1max;u2min=min(u2);u2max=max(u2);u3min=u2max;u3max=max(u3);u4min=u3max;u4max=96.6;x1=[];x2=[];x3=[];x4=[];y=3401;for j=1:3401if u3(j)<=u3minx1=[j];endendq1=max(x1);ua3=u3(q1:y);a3=b3(q1:y);for k=1:3401if u4(k)<=u4minx2=[k];elseif u4(k)<=u4maxx3=[k];endendq2=max(x2);q3=max(x3);ua4=u4(q2:q3);a4=b4(q2:q3);s2=trapz(b2,u2);%二挡运行时间s3=trapz(ua3,a3);%三挡运行时间s4=trapz(ua4,a4);%四挡运行时间s=[s2s3s4];t=sum(s)*1000/3600;%总时间disp('t=');disp(t);gtext(['t='num2str(t)'s']);Pe1=Ft1.*u1./3600;%各档对应转速下的功率Pe2=Ft2.*u2./3600;Pe3=Ft3.*u3./3600;Pe4=Ft4.*u4./3600;endPemax=max([max(Pe1)max(Pe2)max(Pe3)max(Pe4)]);%发动机最大功率disp('90％负荷');Pe=0.9*Pemax%90％负荷rou=0.7;Ua=0.377*r*3401./(i0.*ig4(4));B=[1233.9-84.4782.9788-0.0474490.00028230];%负荷特性曲线拟合公式系数，对应n=3403r/min for i=1:5b=B*[1Pe Pe^2Pe^3Pe^4]';%根据拟合公式计算比油耗Qs4(i)=Pe.*b./(1.02.*Ua(i).*rou.*g);%百公里油耗量enddisp('Qs4=');disp(Qs4);tt=[48.423546.351141.555733.379630.5564];figure(2)plot(Qs4,tt,'*')hold onxi=25:0.001:31;pp=interp1(Qs4,tt,xi,'cubic');plot(xi,pp)title('燃油经济性-加速时间曲线');xlabel('燃油经济性--等速百公里油耗(L/100km)');ylabel('动力性--原地起步加速时间(s)');gtext('5.17')gtext('5.43')gtext('5.83')gtext('6.17')gtext('6.33')%5挡发动机for i=1:5Ft1=Tq*ig5(1)*i0(i)/r;%各档对应转速下的驱动力Ft2=Tq*ig5(2)*i0(i)/r;Ft3=Tq*ig5(3)*i0(i)/r;Ft4=Tq*ig5(4)*i0(i)/r;Ft5=Tq*ig5(5)*i0(i)/r;u1=0.377*r*n/(i0(i)*ig5(1));%各档对应转速下的车速u2=0.377*r*n/(i0(i)*ig5(2));u3=0.377*r*n/(i0(i)*ig5(3));u4=0.377*r*n/(i0(i)*ig5(4));u5=0.377*r*n/(i0(i)*ig5(5));F1=m*g*f+CDA*u1.^2/21.15;%各档对应转速下的行驶阻力F2=m*g*f+CDA*u2.^2/21.15;F3=m*g*f+CDA*u3.^2/21.15;F4=m*g*f+CDA*u4.^2/21.15;F5=m*g*f+CDA*u5.^2/21.15;delta1=1+(Iw+If.*ig5(1).^2.*i0(i).^2.*eta_t)./(m.*r.^2);%汽车旋转质量换算系数delta2=1+(Iw+If.*ig5(2).^2.*i0(i).^2.*eta_t)./(m.*r.^2);delta3=1+(Iw+If.*ig5(3).^2.*i0(i).^2.*eta_t)./(m.*r.^2);delta4=1+(Iw+If.*ig5(4).^2.*i0(i).^2.*eta_t)./(m.*r.^2);delta5=1+(Iw+If.*ig5(5).^2.*i0(i).^2.*eta_t)./(m.*r.^2);a1=(Ft1-F1)/(delta1*m);%加速度a2=(Ft2-F2)/(delta2*m);a3=(Ft3-F3)/(delta3*m);a4=(Ft4-F4)/(delta4*m);a5=(Ft5-F5)/(delta5*m);b1=1./a1;%加速度倒数b2=1./a2;b3=1./a3;b4=1./a4;b5=1./a5;figure(3)subplot(2,3,i)plot(u1,b1,u2,b2,u3,b3,u4,b4,u5,b5);title('加速度倒数-速度曲线');xlabel('u')ylabel('1/a')gtext('1/a1')gtext('1/a2')gtext('1/a3')gtext('1/a4')gtext('1/a5')gtext(['i0='num2str(i0(i))]);%计算采用2挡起步加速至96.6km/h的原地起步加速时间u1min=min(u1);u1max=max(u1);u2min=u1max;u2min=min(u2);u2max=max(u2);u3min=u2max;u3max=max(u3);u4min=u3max;u4max=max(u4);u5min=u4max;u5max=96.6;x1=[];x2=[];x3=[];x4=[];x5=[];y=3401;for j=1:3401if u3(j)<=u3minx1=[j];endendq1=max(x1);ua3=u3(q1:y);a3=b3(q1:y);for k=1:3401if u4(k)<=u4minx2=[k];endendq2=max(x2);ua4=u4(q2:y);a4=b4(q2:y);for l=1:3401if u5(l)<=u5minx3=[l];elseif u5(l)<=u5maxx4=[l];endendq2=max(x2);q3=max(x3);q4=max(x4);ua5=u5(q3:q4);a5=b5(q3:q4);s2=trapz(b2,u2);%二挡运行时间s3=trapz(ua3,a3);%三挡运行时间s4=trapz(ua4,a4);%四挡运行时间s5=trapz(ua5,a5);%五挡运行时间s=[s2s3s4s5];t=sum(s)*1000/3600;%总时间disp('t=');disp(t);gtext(['t='num2str(t)'s']);Pe1=Ft1.*u1./3600;%各档对应转速下的功率Pe2=Ft2.*u2./3600;Pe3=Ft3.*u3./3600;Pe4=Ft4.*u4./3600;Pe5=Ft5.*u5./3600;endPemax=max([max(Pe1)max(Pe2)max(Pe3)max(Pe4)max(Pe4)]);%发动机最大功率disp('90％负荷');Pe=0.9*Pemax%90％负荷rou=0.7;Ua=0.377*r*3401./(i0.*ig4(4));B=[1233.9-84.4782.9788-0.0474490.00028230];%负荷特性曲线拟合公式系数，对应n=3403r/min for i=1:5b=B*[1Pe Pe^2Pe^3Pe^4]';%根据拟合公式计算比油耗Qs5(i)=Pe.*b./(1.02.*Ua(i).*rou.*g);%百公里油耗量enddisp('Qs5=');disp(Qs5);tt=[93.973058.283444.678445.366744.9793];figure(2)plot(Qs5,tt,'*')hold onxi=30:0.001:100;pp=interp1(Qs5,tt,xi,'cubic');plot(xi,pp)%title('燃油经济性-加速时间曲线');%xlabel('燃油经济性--等速百公里油耗(L/100km)');%ylabel('动力性--原地起步加速时间(s)');gtext('5.17')gtext('5.43')gtext('5.83')gtext('6.17')gtext('6.33')运行结果：t=48.4235t=46.3511t=41.5557t=33.3796t=30.556490％负荷Pe=55.5647Qs4=25.103126.365628.307829.958730.7356图1加速度倒数-速度曲线图2燃油经济性-加速时间曲线4.3%空载(no load)—1；满载(full load)—2m1=4080;m2=9290;hg1=0.845;hg2=1.170;L=3.950;a1=2.100;a2=2.950;b1=L-a1;b2=L-a2;beta=0.38;%利用附着系数与制动强度的关系曲线z=0:0.01:1;phi_f1=L*beta.*z./(b1+z*hg1);%前轮利用附着系数(空载)phi_r1=L*(1-beta).*z./(a1-z*hg1);%空载时后轮利用附着系数(空载)phi_f2=L*beta.*z./(b2+z*hg2);%前轮利用附着系数(满载)phi_r2=L*(1-beta).*z./(a2-z*hg2);%后轮利用附着系数(满载)phi=z;%最理想情况figure(1);plot(z,phi_f1,'g--',z,phi_f2,'g-',z,phi_r1,'b--',z,phi_r2,'b-',z,phi,'k-'); axis([0,1,0,1]);grid on;box off;axis square;title('利用附着系数与制动强度的关系曲线');xlabel('制动强度z/g');ylabel('利用附着系数{\phi}');%ECE法规要求界限hold on;z1=0.2:0.01:0.8;z2=0.15:0.01:0.3;z3=0.3:0.01:1;phi1=(z1+0.07)/0.85;phi2=z2-0.08;phi3=z2+0.08;phi4=(z3-0.02)/0.74;plot(z1,phi1,'r-.',z2,phi2,'r-.',z2,phi3,'r-.',z3,phi4,'r-.'); legend('{\phi}_{f}(空载)','{\phi}_{f}(满载)','{\phi}_{r}(空载)',...'{\phi}_{r}(满载)','{\phi}=z','ECE法规要求界限');legend('Location','northwest');%制动效率曲线Ef1=z./phi_f1*100;Er1=z./phi_r1*100;Ef2=z./phi_f2*100;Er2=z./phi_r2*100;figure(2);plot(phi_f1,Ef1,'k',phi_r1,Er1,'k',phi_f2,Ef2,'b',phi_r2,Er2,'b'); axis([0,1,0,100]);grid on;box off;axis square;title('制动效率曲线');xlabel('附着系数{\phi}');ylabel('制动效率(%)');text(0.25,90,'E_f');text(0.6,88,'满载');text(0.8,90,'E_r');text(0.7,65,'空载');text(0.72,72,'E_r');程序运行结果如下：图1图2 5.11L=3.048;K=0.0024;i=20;ua=0:0.01:120;f=1/L.*(1./(3.6./ua+K.*ua./3.6));syms uF(u)=(1/L.*(1./(1./u+K.*u)))./i;plot(ua,f,'b-');xlabel('u_a(km/h)');ylabel('\omega_r/\delta)_s');F(22.35)运行结果：ans=0.1667386094155036.2f=0.1:0.1:100;Gq1=5.12*10^(-5)./f.^2;Gq2=2.02*10^(-3).*f./f;Gq3=7.98*10^(-2).*f.^2;figure(1);loglog(f,Gq1,'b-');box off;axis square;title('位移功率谱密度'); xlabel('f/Hz');ylabel('G_q(f)/(m^2s)');figure(2);loglog(f,Gq2,'b-');box off;axis square;title('速度功率谱密度'); xlabel('f/Hz');ylabel('G_q_''(f)/(m^2/s)');figure(3);loglog(f,Gq3,'b-');box off;axis square;title('加速度功率谱密度'); xlabel('f/Hz');ylabel('G_q_''_''(f)/(m^2/s^3)');运行结果如下：6.4f1=0.1:0.01:1;f2=1:0.01:10;G1=0.010106.*f1.^2;G2=0.010106.*f2./f2;loglog(f1,G1,'b-');hold onloglog(f2,G2,'b-');title('车身加速度的功率谱密度'); xlabel('f/Hz');ylabel('G_z_''_''(f)/m^-^1'); axis([0.1,10,0.0001,0.1]);运行结果如下：6.51.计算说明1）①幅频特性()122222114z q λζλγ⎡⎤-+⎢=⎢⎥∆⎣⎦，其中0ωλω=，()2222222111141λλγζλγλμμ⎡⎤⎡⎤⎛⎫⎛⎫∆=-+--+-⎢⎥⎢⎥ ⎪ ⎪⎝⎭⎝⎭⎣⎦⎣⎦；()12222222211414z z ζλλζλ⎡⎤+⎢=⎢⎥-+⎣⎦；()1222222221414s s s s s p z ζλλζλ⎡⎤+⎢=⎢⎥-+⎣⎦，其中s s ωλω=；②均方根谱()()221~j 2z q H f f ωππ=；()()222~j 2z qH f f ωππ== ；()()2221~1j 2p qz z H f f z qωππ== ；③其他值()0.5360q q G f df σ⎡⎤=⎢⎥⎣⎦⎰ ，()0.536110z z G f df σ⎡⎤=⎢⎥⎣⎦⎰ ，()0.536220z z G f df σ⎡⎤=⎢⎥⎣⎦⎰ ，()0.5360a a G f df σ⎡⎤=⎢⎥⎣⎦⎰()()()()()()13622w 012222412.5362202412.5d 12.5 0.5d d 1d d 4a a a a a a W f G f f f G f f G f f G f f G f f f ⎡⎤=⎢⎥⎣⎦⎡⎤⎛⎫⎛⎫=+++⎢⎥ ⎪ ⎪⎝⎭⎢⎥⎝⎭⎣⎦⎰⎰⎰⎰⎰其中，()()()()()0.5 0.52 2441 412.512.512.536f f f W f f f f <<⎧⎪⎪<<⎪=⎨<<⎪⎪<<⎪⎩()aw w 020lg L a =2）见1）中计算式3）()()()22004q q q G f G G n n uωπ== ()()()1122363622d ~00d d dd f q q f q f H f G f f G f f q σ⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎰⎰ ，其中122d 12f q f γλπ⎡⎤=⎢⎥∆⎣⎦ ()()()1122363622d //~00d d dd F Gq q F G q F H f G f f G f f Gq σ⎡⎤⎡⎤==⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦⎰⎰ ，其中122222d 1412F f Gq g λζλμπγ⎡⎤⎛⎫-+⎢⎥⎪+⎢⎥⎝⎭=⎢⎥∆⎢⎥⎢⎥⎣⎦其余见1）中计算式2.程序清单1）f0=1.5;zeta=0.25;gamma=9;mu=10;fs=3;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega0=2*pi*f0;omega_s=2*pi*fs;lambda=omega./omega0;lambda_s=omega./omega_s;%计算并绘制幅频特性delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;p_z2=((1+(2*zeta_s.*lambda_s).^2)./((1-lambda_s.^2).^2+(2*zeta_s.*...lambda_s).^2)).^0.5;figure(1)loglog(f,z1_q,'b-',0:0.1:100,(0:0.1:100)./(0:0.1:100),'g-');axis([0.1,100,0.1,10]);grid on ;title('z_1~q 的幅频特性');xlabel('激振频率f/Hz');ylabel('$$|\frac{z_1}{q}|$$','Interpreter','latex');gtext('0:1');figure(2)loglog(f,z2_z1,'b-',0:0.1:1,(0:0.1:1)./(0:0.1:1),'g-',1:0.1:100,1./...(1:0.1:100),'g-');axis([0.1,100,0.01,10]);grid on;title('z_2~z_1的幅频特性');xlabel('激振频率f/Hz');ylabel('$$|\frac{z_2}{z_1}|$$','Interpreter','latex'); gtext('0:1');gtext('-1:1');figure(3)loglog(f,p_z2,'b-');axis([0.1,100,0.01,10]);grid on;title('p~z_2的幅频特性');xlabel('激振频率f/Hz');ylabel('$$|\frac{p}{z_2}|$$','Interpreter','latex');%计算并绘制均方根谱sqrt_Gq=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u);sqrt_Gz1=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z1_q;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sqrt_Ga=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*p_z2.*z2_z1.*z1_q; figure(4)loglog(f,sqrt_Gz1,'b-');grid on;title('车轮位移均方根谱');xlabel('激振频率f/Hz');ylabel('$\sqrt{G_{z''''1}(f)}$','Interpreter','latex'); figure(5)loglog(f,sqrt_Gz2,'b-');grid on;title('车身位移均方根谱');xlabel('激振频率f/Hz');ylabel('$\sqrt{G_{z''''2}(f)}$','Interpreter','latex'); figure(6)loglog(f,sqrt_Ga,'b-');grid on;title('传至人体的位移均方根谱');xlabel('激振频率f/Hz');ylabel('$\sqrt{G_a(f)}$','Interpreter','latex');%计算其它值sigma_q=sqrt(trapz(f,sqrt_Gq.^2));%路面不平度加速度均方根值sigma_z1=sqrt(trapz(f,sqrt_Gz1.^2));%车轮加速度均方根值sigma_z2=sqrt(trapz(f,sqrt_Gz2.^2));%车身加速度均方根值sigma_a=sqrt(trapz(f,sqrt_Ga.^2));%传至人体的加速度均方根值for i=1:Nif f(i)<=2W(i)=0.5;elseif f(i)<=4W(i)=f(i)/4;elseif f(i)<=12.5W(i)=1;elseW(i)=12.5/f(i);endendaw=sqrt(trapz(f,W.^2.*sqrt_Ga.^2));%加权加速度均方根值a0=10^(-6);Law=20*log10(aw/a0);%加权振级format shortdisp('路面不平度加速度均方根值=');disp(sigma_q);disp('车轮加速度均方根值=');disp(sigma_z1);disp('车身加速度均方根值=');disp(sigma_z2);disp('传至人体的加速度均方根值=');disp(sigma_a);disp('加权加速度均方根值=');disp(aw);disp('加权振级=');disp(Law);2）%随fs变化f0=1.5;zeta=0.25;gamma=9;mu=10;fs=1.5:0.01:6;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega0=2*pi*f0;lambda=omega./omega0;for i=1:length(fs)omega_s=2*pi.*fs(i);lambda_s=omega./omega_s;delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;p_z2=((1+(2*zeta_s.*lambda_s).^2)./((1-lambda_s.^2).^2+(2*zeta_s.*...lambda_s).^2)).^0.5;sqrt_Gq=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u);sqrt_Gz1=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z1_q;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sqrt_Ga=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*p_z2.*z2_z1.*z1_q;for j=1:Nif f(j)<=2W(j)=0.5;elseif f(j)<=4W(j)=f(j)/4;elseif f(j)<=12.5W(j)=1;elseW(j)=12.5/f(j);endendaw(i)=sqrt(trapz(f,W.^2.*sqrt_Ga.^2));enda0=10^(-6);Law=20*log10(aw/a0);figure(1)plot(fs,aw);grid ontitle('a_w随f_s的变化曲线');xlabel('f_s/Hz');ylabel('a_w/m·s^-^2');figure(2)plot(fs,Law);grid ontitle('L_aw随f_s的变化曲线');xlabel('f_s/Hz');ylabel('L_a_w/dB');%随zeta_s变化f0=1.5;zeta=0.25;gamma=9;mu=10;fs=3;zeta_ss=0.125:0.001:0.5;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega0=2*pi*f0;omega_s=2*pi.*fs;lambda=omega./omega0;lambda_s=omega./omega_s;for i=1:length(zeta_ss)zeta_s=zeta_ss(i);delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;p_z2=((1+(2*zeta_s.*lambda_s).^2)./((1-lambda_s.^2).^2+(2*zeta_s.*...lambda_s).^2)).^0.5;sqrt_Gq=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u);sqrt_Gz1=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z1_q;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sqrt_Ga=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*p_z2.*z2_z1.*z1_q;for j=1:Nif f(j)<=2W(j)=0.5;elseif f(j)<=4W(j)=f(j)/4;elseif f(j)<=12.5W(j)=1;elseW(j)=12.5/f(j);endendaw(i)=sqrt(trapz(f,W.^2.*sqrt_Ga.^2));enda0=10^(-6);Law=20*log10(aw/a0);figure(1)plot(zeta_ss,aw);grid ontitle('a_w随\zeta_s的变化曲线');xlabel('\zeta_s');ylabel('a_w/m·s^-^2');figure(2)plot(zeta_ss,Law);grid ontitle('L_aw随\zeta_s的变化曲线');xlabel('\zeta_s');ylabel('L_a_w/dB');3）%随f0变化f0=0.25:0.01:3;zeta=0.25;gamma=9;mu=10;fs=3;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega_s=2*pi*fs;lambda_s=omega./omega_s;Gq_f=4*pi^2*Gq_n0*n0^2*u;g=9.8;for i=1:length(f0)omega0=2*pi*f0(i);lambda=omega./omega0;delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sigma_z2(i)=sqrt(trapz(f,sqrt_Gz2.^2));fd_q=gamma*lambda.^2./(2*pi.*f).*(1./delta).^0.5;sigma_fd(i)=sqrt(trapz(f,Gq_f.*fd_q.^2));Fd_Gq=2*pi.*f.*gamma./g.*(((lambda.^2./(1+mu)-1).^2+4*zeta.^2*lambda.^2)..../delta).^0.5;sigma_FdG(i)=sqrt(trapz(f,Gq_f.*Fd_Gq.^2));endm=find(f0==1.5);sgm_z2=sigma_z2(m);sgm_fd=sigma_fd(m);sgm_FdG=sigma_FdG(m);sm_z2=20*log10(sigma_z2/sgm_z2);sm_fd=20*log10(sigma_fd/sgm_fd);sm_FdG=20*log10(sigma_FdG/sgm_FdG);plot(f0,sm_z2,f0,sm_fd,f0,sm_FdG);axis([0.25,3,-25,10]);grid on;title('各响应量均方根值随f_0变化的曲线');xlabel('f_0/Hz');ylabel('\sigma_z_''_''_2,\sigma_f_d,\sigma_F_d_/_G/dB'); legend('\sigma_z_''_''_2','\sigma_f_d','\sigma_F_d_/_G');legend('Location','southeast');%随zeta变化f0=1.5;zeta0=0.125:0.001:0.5;gamma=9;mu=10;fs=3;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega_s=2*pi*fs;omega0=2*pi*f0;lambda_s=omega./omega_s;lambda=omega./omega0;Gq_f=4*pi^2*Gq_n0*n0^2*u;g=9.8;for i=1:length(zeta0)zeta=zeta0(i);delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sigma_z2(i)=sqrt(trapz(f,sqrt_Gz2.^2));fd_q=gamma*lambda.^2./(2*pi.*f).*(1./delta).^0.5;sigma_fd(i)=sqrt(trapz(f,Gq_f.*fd_q.^2));Fd_Gq=2*pi.*f.*gamma./g.*(((lambda.^2./(1+mu)-1).^2+4*zeta.^2*lambda.^2)..../delta).^0.5;sigma_FdG(i)=sqrt(trapz(f,Gq_f.*Fd_Gq.^2));endm=find(zeta0==0.25);sgm_z2=sigma_z2(m);sgm_fd=sigma_fd(m);sgm_FdG=sigma_FdG(m);sm_z2=20*log10(sigma_z2/sgm_z2);sm_fd=20*log10(sigma_fd/sgm_fd);sm_FdG=20*log10(sigma_FdG/sgm_FdG);plot(zeta0,sm_z2,zeta0,sm_fd,zeta0,sm_FdG);axis([0.125,0.5,-4,4]);grid on;title('各响应量均方根值随\zeta变化的曲线');xlabel('\zeta');ylabel('\sigma_z_''_''_2,\sigma_f_d,\sigma_F_d_/_G/dB'); legend('\sigma_z_''_''_2','\sigma_f_d','\sigma_F_d_/_G');%随gamma变化f0=1.5;zeta=0.25;gamma0=4.5:0.1:18;mu=10;fs=3;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega_s=2*pi*fs;omega0=2*pi*f0;lambda_s=omega./omega_s;lambda=omega./omega0;Gq_f=4*pi^2*Gq_n0*n0^2*u;g=9.8;for i=1:length(gamma0)gamma=gamma0(i);delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sigma_z2(i)=sqrt(trapz(f,sqrt_Gz2.^2));fd_q=gamma*lambda.^2./(2*pi.*f).*(1./delta).^0.5;sigma_fd(i)=sqrt(trapz(f,Gq_f.*fd_q.^2));Fd_Gq=2*pi.*f.*gamma./g.*(((lambda.^2./(1+mu)-1).^2+4*zeta.^2*lambda.^2)..../delta).^0.5;sigma_FdG(i)=sqrt(trapz(f,Gq_f.*Fd_Gq.^2));endm=find(gamma0==9);sgm_z2=sigma_z2(m);sgm_fd=sigma_fd(m);sgm_FdG=sigma_FdG(m);sm_z2=20*log10(sigma_z2/sgm_z2);sm_fd=20*log10(sigma_fd/sgm_fd);sm_FdG=20*log10(sigma_FdG/sgm_FdG);plot(gamma0,sm_z2,gamma0,sm_fd,gamma0,sm_FdG);axis([4.5,18,-6,6]);grid on;title('各响应量均方根值随\gamma变化的曲线');xlabel('\gamma');ylabel('\sigma_z_''_''_2,\sigma_f_d,\sigma_F_d_/_G/dB'); legend('\sigma_z_''_''_2','\sigma_f_d','\sigma_F_d_/_G');legend('Location','southeast');%随mu变化f0=1.5;zeta=0.25;gamma=9;mu0=5:0.1:20;fs=3;zeta_s=0.25;u=20;Gq_n0=2.56*10^(-8);n0=0.1;delta_f=0.2;N=180;f=delta_f*[1:N];omega=2*pi*f;omega_s=2*pi*fs;omega0=2*pi*f0;lambda_s=omega./omega_s;lambda=omega./omega0;Gq_f=4*pi^2*Gq_n0*n0^2*u;g=9.8;for i=1:length(mu0)mu=mu0(i);delta=((1-lambda.^2).*(1+gamma-1./mu.*lambda.^2)-1).^2+4*zeta^2.*lambda.^...2.*(gamma-(1./mu+1).*lambda.^2).^2;z1_q=gamma.*(((1-lambda.^2).^2+4*zeta.^2.*lambda.^2)./delta).^0.5;z2_z1=((1+4*zeta.^2.*lambda.^2)./((1-lambda.^2).^2+4*zeta.^2.*lambda.^2))....^0.5;sqrt_Gz2=4*pi^2.*f.*sqrt(Gq_n0*n0^2*u).*z2_z1.*z1_q;sigma_z2(i)=sqrt(trapz(f,sqrt_Gz2.^2));fd_q=gamma*lambda.^2./(2*pi.*f).*(1./delta).^0.5;sigma_fd(i)=sqrt(trapz(f,Gq_f.*fd_q.^2));Fd_Gq=2*pi.*f.*gamma./g.*(((lambda.^2./(1+mu)-1).^2+4*zeta.^2*lambda.^2)..../delta).^0.5;sigma_FdG(i)=sqrt(trapz(f,Gq_f.*Fd_Gq.^2));endm=find(mu0==10);sgm_z2=sigma_z2(m);sgm_fd=sigma_fd(m);sgm_FdG=sigma_FdG(m);sm_z2=20*log10(sigma_z2/sgm_z2);sm_fd=20*log10(sigma_fd/sgm_fd);sm_FdG=20*log10(sigma_FdG/sgm_FdG);plot(mu0,sm_z2,mu0,sm_fd,mu0,sm_FdG);axis([5,20,-2,2]);grid on;title('各响应量均方根值随\mu变化的曲线');xlabel('\mu');ylabel('\sigma_z_''_''_2,\sigma_f_d,\sigma_F_d_/_G/dB');legend('\sigma_z_''_''_2','\sigma_f_d','\sigma_F_d_/_G');3.结果分析1）路面不平度加速度均方根值=0.3523车轮加速度均方根值=0.2391车身加速度均方根值=0.0168传至人体的加速度均方根值=0.0161加权加速度均方根值=0.0100加权振级=80.0287分析：根据课本中表6-2，a w=0.0100<0.315，L aw=80.0287<110，故乘客没有不舒适。