MATLAB软件与基础数学实验

MATLAB软件与基础数学实验 实验1 MATLAB基本特性与基本运算 例1-1 求[12+2×(7-4)]÷32的算术运算结果。 ► （12+2*（7-4））/3^2 ◄ ans = 2 例1-2 计算5!，并把运算结果赋给变量y ► y=5*4*3*2*1 ◄ y =120 例1-3 ► sqrt(2) % 计算2开平方 ◄ ans = 1.4142 例1-4 ► x = sqrt(2)； % 计算2开平方并赋值给变量x（不显示） ► x % 查看x的赋值情况 ◄ x = 1.4142 例1-5 设 ，计算 的值。 ► a=pi/180*(-24); %转换为弧度值且不显示 ► b=pi/180*75; %转换为弧度值且不显示 ► z=sin(abs(a)+abs(b))/sqrt(tan(abs(a+b))) %计算结果并显示 ◄ z =0.8888 例1-6 设三角形三边长为 ，求此三角形的面积。 ► a=4;b=3;c=2; %输入边长值且不显示 ► s=(a+b+c)/2; ► A=s*(s-a)*(s-b)*(s-c); %计算面积平方且不显示 ► A=sqrt(A) %计算面积并显示 ◄ A = 2.9047 例1-7 设 ， ，计算 ， 。 ► A=[1,2,3;4,5,6;1,0,1]; ► B=[-1 2 0;1 1 3;2 1 1]; ► ◄ans= 0 4 3 5 6 9 3 1 2 ► ◄Ans=7 7 9 13 19 21 1 3 1 ► % det为求方阵的行列式命令 ◄ans = -6 ► %inv为方阵的求逆命令 ◄ans = -0.8333 0.3333 0.5000 -0.3333 0.3333 -1.0000 0.8333 -0.3333 0.5000 例1-8 显示上例中矩阵A的第2行第3列元素，并对其进行修改. ► A(2,3) ◄ A(2,3) = 6 若想把该元素改为-1，只要输入下列语句： ► A(2,3)= -1； 例1-9 分别画出函数 和 在区间[-6 ,6 ]上的图形。 ► x= (-6 : 0.1 : 6)*pi; % 从-6pi到6pi以0.1pi为步长生成向量x ► y=x.^2 .* cos(x); % 产生与x对应的函数值向量y(两向量对应元素乘积，用.*) ► z=sin(x) ./ (x+eps); % 产生与x对应的函数值向量z(两向量对应元素相除，用./) ► subplot(1,2,1) % 分图形窗口为1行2列，并在第一个子窗中绘图 ► plot(x,y,'linewidth',2) % 画函数y的曲线,默认为蓝色（参看实验2） ► grid %在第一个子窗中加坐标网格 ► subplot(1,2,2) %在第二个子窗中绘图 ► plot(x,z,'linewidth',2) % 画函数z的曲线,默认为蓝色（参看实验2） ► grid %在第二个子窗中加坐标网格 例1-10 试求方程组 的解。 ► a=[1,2,1;4,2,-6;-1,0,2]； % 输入系数矩阵a ► b=[2;3;4]； % 输入右端列向量b ► d=det(a) %求系数矩阵的行列式 ◄ d= 2 ► c=inv(a) %求系数矩阵的逆阵 ◄ c = 2.0000 -2.0000 -7.0000 -1.0000 1.5000 5.0000 1.0000 -1.0000 -3.0000 ► x=c*b %矩阵左逆乘，结果为方程组的解 ◄ x = -30.0000 22.5000 -13.0000 ► X=a\b %用\除法直接求方程组的解X（与上述x相同） ◄ X = -30.0000 22.5000 -13.0000 ► disp([a,b,x]) %显示增广矩阵及解向量 ◄ 1.0000 2.0000 1.0000 2.0000 -30.0000 4.0000 2.0000 -6.0000 3.0000 22.5000 -1.0000 0 2.0000 4.0000 -13.0000 例1-11 试求矩阵方程 的解。 ► a=[1,2,1;4,2,-6;-1,0,2] ; % 输入系数矩阵a ► b=[1 2 3;1 1 1] ; % 输入右端矩阵b ► X=b/a % 用/除法直接求方程组的解X ◄ X = 3.0000 -2.0000 -6.0000 2.0000 -1.5000 -5.0000 例1-12 建立同时计算 ， 的函数。即任给a,b,n三个数，返回y1,y2. function [y1 , y2]=fun1(a , b , n) % fun1 is a function used by DEMO y1=(a+b)^n, y2=(a-b)^n % Copyright by XJTU y1=(a+b).^n ; y2=(a-b).^n; 例1-13 设 ，试画出在[0,2]上的曲线段。 ► x=0 : 0.01 : 2; %生成自变量x ► y=1 ./ ((x-0.3) .^2+0.01)+1 ./ ((x-0.9) .^2+0.04)-6; %生成函数值y，注意点运算 ► plot(x,y,'linewidth',2) %画函数曲线 ► grid %加坐标网格 ► f=inline(' 1 ./ ((x-0.3) .^2+0.01)+1 ./ ((x-0.9) .^2+0.04)-6 '); %生成数值函数f(x) ► fplot(f,[0,2]) % 画函数f在[0,2]上的曲线 ► grid % 加坐标网格 例如：对于例题1-13中所定义的f(x)，求其零点c. ► f=inline(' 1 ./ ((x-0.3) .^2+0.01)+1 ./ ((x-0.9) .^2+0.04)-6 '); %生成数值函数f(x) ► c=fzero(f , [0,2]) % 求函数f在[0,2]上的零点c，此处要求f(0)f(2)<0 ◄ c= 1.2995 ► fzero(f , 1) % 求函数f在x=1附近的零点 ◄ ans = 1.2995 例如：求一元函数最小值（fminbnd命令） ► fy=inline('1./((x-0.3).^2+0.01)+1./((x-0.9).^2+0.04)-6'); ► [xmin,fmim]=fminbnd(fy,0.2,0.8) %函数fy在[0.2,0.8]上最小值点及最小值 ◄ xmin = 0.6370 ◄ fmim = 11.2528 ► ff=inline('-1./((x-0.3).^2+0.01)-1./((x-0.9).^2+0.04)+6'); %函数ff=-fy ► [x,y]=fminbnd(ff,0.2,0.8); %函数ff在[0.2,0.8]上最小值点及最小值 ► xmax=x ◄ xmax = 0.3004 ► fmax=-y ◄ fmax = 96.5014 例如：求例题1-13中所定义f(x)在[0,1]上的定积分 . ► f=inline('1./((x-0.3).^2+0.01)+1./((x-0.9).^2+0.04)-6'); ► I=quad(f,0,1) % 求f(x)在[0,1]上定积分 ◄ I = 29.8583 例1-14 求二重积分 及三重积分 。 ► g=inline('x.*y', 'x','y'); % 建立二元函数g(x,y)=xy ► I=dblquad(g,0,1,1,2) % 求g(x,y)在[0,1] ×[1,2]上的二重积分 ◄ I = 0.7500 ► h=inline('x.*exp(y)+z.^2', 'x','y','z'); % 建立三元函数 ► I=triplequad(h,0,1,0,1,0,1) % 求h(x,y,z)在[0,1] ×[0,1] ×[0,1]上的三重积分 ◄ I = 1.1925 例1-15 已知 ，设该曲线在区间[0,x]上所围曲边梯形面积为s，试求当s分别为5，10时的x的值。 分. (1) 对于s=5 ► f=inline('1/4*x^4-5/3*x^3+3*x^2+5*x-5') ; %建立函数 ► x=fzero(f,[0,5]) %求解方程 在[0,5]上的根 ◄ x =0.7762 (2) 对于s=10 ►g=inline('1/4*x^4-5/3*x^3+3*x^2+5*x-10'); ► x=fzero(g,[0,10]) %求解方程 在[0,10]上的根 ◄ x =1.5179 例1-16 利用MATLAB命令求解无理数的近似值。 (1) 用函数零点命令（fzero）求无理数 的近似值； (2) 用定积分计算命令（trapz,quad,quadl）求无理数 的近似值。 （提示：e =2.7182818284…， =0.6931471806…） (1) 无理数 可以看成是方程 在x=2附近的实根，于是可以用fzero来求解。 ► f=inline('log(x)-1') ; %建立函数 ► x0=fzero(f,2); %求解方程 在x=2附近的根 ► e =vpa(x0,10) %显示x0小数点后10位 ◄ e =2.7182818284 (2) 由于无理数 ，于是可以用trapz,quad,quadl命令分别来求解。 用梯形法（trapz）近似计算 ► X=0:0.01:1; %产生[0,1]区间上的划分向量 ► Y=1./(1+X); %求对应的分点处的函数值向量 ► a=trapz(X,Y); %求用梯形法求出积分近似值 ► ln2 =vpa(a,10) %显示a小数点后10位 ◄ ln2 =0.6931534305 （注意：已精确到小数点后4位） 用高阶方法（quad,quadl）近似计算 ► f=inline('1./(1+x)'); %建立被积函数f(x) ► a=quad(f,0,1); %用辛浦生方法求f在[0,1]上的积分近似值 ► ln2 =vpa(a,10) %显示a小数点后10位 ◄ ln2 =0.6931471999 （注意：已精确到小数点后7位） ► a=quadl(f,0,1); %用高阶方法求f在[0,1]上的积分近似值 ► ln2 =vpa(a,10) %显示a小数点后10位 ◄ ln2 =0.6931471861 （注意：已精确到小数点后9位） 例1-17 求极限 。 ►syms h ► fx= sym ('(sin(x+h)-sin(x))/h') ; %建立符号函数fx ► limit(fx,h,0) %求fx : h->0的极限 ◄ ans=cos(x) 例1-18：设 ，求 ► syms x y n %声明符号变量，注意变量间必须用空格分开 ► fx=x^n*y+sin(y); %建立符号函数 ► diff(fx) %对变量x(默认)求一阶导数（偏导数） ◄ ans =x^n*n/x*y 即 ► diff(fx, y) %对变量y求一阶导数（偏导数） ◄ ans =x^n+cos(y) ► diff(fx, y, 2) %对变量y求二阶导数（偏导数） ◄ ans =-sin(y) ► diff(diff(fx,x), y) %先对x求导再对y求导（二阶混合偏导数） ◄ ans = x^n*n/x 即 例1-19：求 ， ， ， ► syms x y z %声明符号变量，注意变量间必须用空格分开 ► f1=x*y/(1+x^2) ; %建立符号函数 ► f2=x+y+z; ► int(f1) %对f1关于变量x(默认)求不定积分 ◄ ans =1/2*y*log(1+x^2) %即 ►syms t ► int(f1,0, t ) %对f1关于变量x(默认)在[0,t]上求定积分 ◄ ans =1/2*log(1+t^2)*y %即 ► int(int(f1,y,0, sqrt(x)),x,0,1 ) %对f1先求对y的积分再求对x的积分（二重积分） ◄ ans =1/2-1/8*pi %即 ► int(int(int(f2,z,0, 1-x-y),y,0,1-x),x,0,1 ) %对f2先对zy的积分再求对x的积分（二重积分） ◄ ans =1/8 级数求和（symsum） ► syms a k ► symsum(1/k,1,inf) %求级数 (ans=inf 即 ) ► symsum(1/(k*(k+1)),1,inf) %求级数 (ans=1) ► symsum(a*1/3^k,k,0,inf) %求级数 (ans= 3/2*a) 泰勒展开（taylor） ► syms x ► fy=1/(1+x+x^2) ► f=taylor(fy) %求fx对自变量x(默认)在x=0点(默认)泰勒展开前6项(默认) ► f=taylor(fy,8,1) %求fx对自变量x(默认)在x=1点泰勒展开式前8项 方程求根（solve） ► fx=sym('a*x^2+b*x+c') ; %建立符号函数 ► solve(fx) %求方程fx=0的符号解 ◄ ans = [ 1/2/a*(-b+(b^2-4*a*c)^(1/2))] [ 1/2/a*(-b-(b^2-4*a*c)^(1/2))] ► syms b ► solve(fx, b ) %求方程fx=0关于变量b的符号解 ◄ ans = -(a*x^2+c)/x 微分方程(组)求解（dsolve） ► dsolve('Dy=5') %求方程y'=5的通解，默认自变量为t ◄ ans = 5*t+C1 ► dsolve('Dy=x', 'x') %求方程y'=x的通解，指定自变量为x ◄ ans =1/2*x^2+C1 ► dsolve('D2y=1+Dy', 'y(0)=1', 'Dy(0)=0') %求方程y''=1+y'满足y(0)=1,y'(0)=0的特解 ◄ ans = -t+exp(t) 即 ► [x,y]=dsolve('Dx=x+y,Dy=2*x') %求方程组 的通解，默认自变量为t ◄ x =1/3*C1*exp(-t)+2/3*C1*exp(2*t)+1/3*C2*exp(2*t)-1/3*C2*exp(-t) y =2/3*C1*exp(2*t)-2/3*C1*exp(-t)+2/3*C2*exp(-t)+1/3*C2*exp(2*t) 即 实验2 MATLAB绘制二维、三维图形 例2-1 在子图形窗口中画出 上正弦、余弦曲线。 ► x=0:0.1*pi:2*pi; %按步长赋值生成x向量 ► y=sin(x); z=cos(x); %生成正弦、余弦函数值y、z向量 ► subplot(2,1,1) %分图形窗口为2行1列，并在第一个子窗中绘图 ► plot(x,y,x,z) %在第一个子窗中画出正弦、余弦曲线 ► subplot(2,1,2) %在第二个子窗中绘图 ► plot(x,y,'k:',x,z,'r-') %在第二个子窗中用不同颜色画两条曲线 ► hold on %保持第二个子窗中绘图 ► plot(x,y,'bo',x,z,'k+') %用'o'和'+'标记曲线上分点 ► hold off %取消图形保持 例2-2 画出 上正弦、余弦曲线并对线型加粗、点型加大，重新定置坐标系以及加注相关说明和注释。 ► x=0:0.1*pi:2*pi; %按步长赋值生成x向量 ► y=sin(x); %生成正弦、余弦函数值y、z向量 ► z=cos(x); ► plot(x,y, 'b-', x,z, 'k .-', ' linewidth ',3, ' markersize ',15) ► axis([-0.2*pi 2.2*pi –1.2 1.2]) %重新设置图形窗口坐标轴范围 ► grid %加注坐标网格 ► xlabel('Variable \it{x}') %标记横坐标轴, \it{x} 表 关于同志近三年现实表现材料材料类招标技术评分表图表与交易pdf视力表打印pdf用图表说话 pdf 示x为斜体 ► ylabel('Variable \it{y}') %标记纵坐标轴 ► title('Sine and Cosine Cruves') %标记图名 ► text(2.5,0.7,'Sin(x)') %在(2.5,0.7)位置，标记曲线名称 ► text(1.5,0.1,'Cos(x)') %在(1.5,0.1)位置，标记曲线名称 ► hold on %图形保持,在同一图形窗口中叠加图形 ► plot([0,2*pi],[0,0], 'r-.') %叠加一条红色的点划直线:(0,0)到(2pi,0) ► hold off %图形保持取消，再画图时将另辟窗口 例2-3 分别在两个图形窗口画出填充一正方形和极坐标方程 的图形。 ► h1=figure; %打开第一个图形窗口，返回其图标识号（句柄）h1 ► x=[0 1 1 0 0]; %闭合图形的顶点横坐标向量 ► y=[0 0 1 1 0]; %闭合图形的顶点纵坐标向量 ► fill(x,y,'y') %填充闭合图形（用黄颜色） ► axis([-1 2 -1 2]) %重新设置坐标轴 ► h2=figure; %打开第二个图形窗口，返回其图标识号（句柄）h2 ► theta=linspace(0,2*pi)； %对theta角的范围进行划分，生成分点向量 ► rho=sin(2*theta).*cos(2*theta); %生成相应极坐标方程的极径rho向量 ► polar(theta,rho,'r') %绘制相应的极坐标方程图形（用红颜色） ► title('Polar plot of sin(2*theta)cos(2*theta)') %添加图形标题 ► set(h2,'linewidth',3) %对第二个窗口中曲线加粗 例2-4在[-2.5,2.5]上画出函数 的直方图和阶梯图。 ► x=linspace(-2.5,2.5,20); %产生横坐标x向量 ► y=exp(-x.*x); %生成函数值向量 ► h1=subplot(1,2,1); %分图形窗口并在第一个子窗中绘图，返回其句柄h1 ► bar(x,y) %画出直方图 ► title(' Bar Chart of a Bell Curve ') %添加图形标题 ► h2= subplot(1,2,2); %在第二个子窗中绘图，返回其句柄h2 ► stairs(x,y) %画出阶梯图 ► title(' Stairs Plot of a Bell Curve ') %添加图形标题 例2-5 采用不同形式（直角坐标、 参数 转速和进给参数表a氧化沟运行参数高温蒸汽处理医疗废物pid参数自整定算法口腔医院集中消毒供应 、极坐标），画出单位圆 的图形。 （1）直角坐标系 ► x=-1:0.01:1； %对x的范围进行划分，生成分点向量 ►y1=sqrt(1-x.^2); %生成上半单位圆的函数值向量 ►y2=-y1; %生成下半单位圆的函数值向量 ►plot(x,y1,x,y2); %同时画出上半圆和下半圆 ►axis equal %让坐标系中两个坐标轴取值相同 （2）参数方程 ► t=0:0.01*pi:2*pi； %对t的范围进行划分，生成分点向量 ►x=cos(t); y=sin(t); %生成单位圆上的函数值向量 ►plot(x,y); %画出单位圆 ►axis equal %让坐标系中两个坐标轴取值相同 （3）极坐标系 ► t=0:0.01*pi:2*pi； %对t的范围进行划分，生成分点向量 ► r=1+0*t; %生成单位圆的极径r向量 ► polar(t,r) %绘制相应的极坐标方程图形 例2-6 画出螺旋线：x=sin(t),y=cos(t),z=t, 上一段曲线。 ► t=0:pi/50:10*pi; %生成参数t数组 ► X=sin(t); %生成螺旋线X数组 ► Y=cos(t); %生成螺旋线Y数组 ► Z=t; %生成螺旋线Z数组 ► plot3(X,Y,Z, 'k-', 'linewidth',3) %画螺旋线 ► grid 例2-7 画出矩形域[-1,1] ×[-1,1]上旋转抛物面： 。 ► x=linspace(-1,1,100); %分割[-1,1]区间生成x ► y=x; %y与x相同 ► [X,Y]=meshgrid(x,y); %生成矩形域[-1,1] ×[-1,1]网格节点坐标矩阵 ► Z=X.^2+Y.^2; %生成 函数值矩阵 ► subplot(1,2,1) ► mesh(X,Y,Z) ; %在第一个子图中画 网格曲面 ► subplot(1,2,2) ► surf(X,Y,Z) ; %在第二个子图中画 光滑曲面 ► shading flat ; %对曲面 平滑并除去网格 例2-8 在圆形域 上绘制旋转抛物面： 。 ► x=linspace(-1,1,300); %分割[-1,1]区间生成x ► y=x; %生成y ► [X,Y]=meshgrid(x,y); %生成矩形域[-1,1]X[-1,1]网格节点坐标矩阵 ► Z=X.^2+Y.^2; %生成 函数值矩阵 ► i=find(Z>1); %找出圆域 之外的函数值（z>1）坐标点i ► Z(i)=NaN; %对圆域 之外的坐标点i处函数值进行“赋空” ► subplot(1,2,1) ► mesh(X,Y,Z) ; %在第一个子图中画 网格曲面 ► subplot(1,2,2) ► surf(X,Y,Z) ; %在第二个子图中画 光滑曲面 ► shading flat ; %对曲面 平滑并除去网格 例2-9 画出 在 上的图形。 ► x=-7.5:0.5:7.5; ► y=x; ► [X,Y]=meshgrid(x,y); ► u=sqrt(X.^2+Y.^2)+eps; %加eps使得u不等于0，保证z有意义 ► Z=sin(u)./u; ► surf(X,Y,Z) 例2-10 有一组实验数据如下表所示，试绘图表示。 时 间 1 2 3 4 5 6 7 8 9 数据1 12.51 13.54 15.60 15.92 20.64 24.53 30.24 50.00 36.34 数据2 9.87 20.54 32.21 40.50 48.31 64.51 72.32 85.98 89.77 数据3 10.11 8.14 14.17 10.14 40.50 39.45 60.11 70.13 40.90 ► t=1:9; ► d1=[12.51 13.54 15.60 15.92 20.64 24.53 30.24 50.00 36.34]; ► d2=[ 9.87 20.54 32.21 40.50 48.31 64.51 72.32 85.98 89.77]; ► d3=[10.11 8.14 14.17 10.14 40.50 39.45 60.11 70.13 40.90]; ► plot(t,d1,'r+-',t,d2,'kx:',t,d3,'b*-','linewidth',2,'markersize',8); ► title('time & data'); ► xlabel('time');ylabel('data'); ► axis([0 10 0 100]); ► text(6.5,25.5,'\leftarrowdata1'); % ' \leftarrow '表示画一左箭头←,且在标识前 ► text(3,43.8,'data2\rightarrow'); % ' \rightarrow '表示画一右箭头→,且在标识后 ► text(4.8,30.5,'\leftarrowdata3'); ► grid 实验3 MATLAB编程介绍与循环结构 例3-1：求n（n=100）个奇数的和：s=1+3+5+…+(2n-1). clear;clc; %清除内存变量，清理命令窗口 n=100; %赋值给定奇数的个数 s=0; %设定存放和的变量s并赋初值0 for i=1:n %定义循环变量i从1到n，以1为步长，即为奇数序号 s=s+(2*i-1); %先计算右端奇数并累加后再赋给左端的变量s fprintf('i=%.0f, s=%.0f\n',i,s) %逐行显示出累加求和的过程 end %循环结构结束 例3-2：求正整数n的阶乘：p=1×2 × 3 × … × n = n!，并求出n=20时的结果。 clear;clc; %清除内存变量，清理命令窗口 n=20; %赋值给定正整数 p=1; %设定存放阶乘的变量p并赋初值1 for i=1:n %定义循环变量i从1到n，以1为步长，即连续正整数 p=p*i; %先计算右端乘积后再赋给左端的变量p fprintf('i=%.0f, p=%.0f\n',i,p) %逐行显示出i! end %循环结构结束 例3-3：根据麦克劳林公式可以得到e≈1+1+1/2!+1/3!+…+1/n!，试求e的近似值。 clear;clc; %清除内存变量，清理命令窗口 n=10; %赋值给定正整数 p=1; %设定存放阶乘的变量p并赋初值1 s=1; %设定存放累加和的变量s并赋初值1 for i=1:n %定义循环变量i从1到n，以1为步长 p=p*i; %先计算右端乘积后再赋给左端的变量p，此时p为i的阶乘 s=s+1/p; %先计算右端阶乘倒数的累加后再赋给左端的变量s fprintf('i=%.0f, s=%.8f\n',i,s) %逐行显示出第i次e的近似值 end %循环结构结束 例3-4：对于数列 ，求其前n项和不超过1000时的n的值及和. clear;clc; %清除内存变量，清理命令窗口 n=0; %设定正整数并赋初值0 s=0; %设定存放累加和的变量s并赋初值0 while s<=1000 %用累加和s与1000进行比较作为循环条件 n=n+1; %改变n为连续正整数 s=s+sqrt(n); %先计算右端开方数的累加后再赋给左端的变量s fprintf('n=%.0f, s=%.4f\n',n,s) %逐行显示正整数及部分和 end %循环结构结束 例3-5：根据e≈1+1+1/2!+1/3!+…+1/n! 求e的近似值，要求精确到 。 clear;clc; %清除内存变量，清理命令窗口 p=1; %设定存放阶乘的变量p并赋初值1 s=1; %设定存放累加和的变量s并赋初值1 r=1; %设定前后两次近似值的误差r并赋初值1 k=0; %设定构造连续正整数的变量k赋初值0又为循环次数 while r>=1.0e-8 %当近似值的精度r没达到 时继续循环 k=k+1; %累计循环次数并作为下一个正整数k p=p*k; %计算k的阶乘p r=1/p; %计算前后两次近似值的误差r s=s+r; %计算e的近似值s fprintf('k=%.0f, s=%.10f\n',k,s) %逐行显示出第k次e的近似值s end %循环结构结束 实验4 MATLAB选择结构与应用实验 例4-1：求任意有限数组a=[a(1),a(2),…,a(n)] 中数值最大的元素M以及所在位置k. function [M,k]=findM(a) %定义函数findM，输入数组a，返回最大元素M及位置k n=length(a); %获取数组的长度即元素的个数n M=a(1); k=1; %将第一个元素作为最大值赋值给M，位置为1; for i=2:n %从第二个元素到最后一个元素依次进行 if a(i)>M %比较后续元素与目前最大值M的大小 M=a(i); k=i; %将数值较大的元素赋值给M，同时保留位置i end %选择结构结束 end %循环结构结束 ► a=[1,2.2,pi,-0.8,3.2,0]; %任意给定一数组 ► [M,k]=findM(a) %调用函数findM ◄ M =3.200 ◄ k =5 例4-2：编写一个函数将百分制成绩转换为优(A)，良(B)，中(C)，差(D)四等级. function jb=dengji(fs) %定义函数dengji，输入分数fs，返回等级A,B,C,D if fs>=90 %判断分数fs是否处在优秀级别上 jb=' A '; %定义为A级 elseif fs>=78 %判断分数fs是否处在良好级别上 jb=' B '; %定义为B级 elseif fs>=60 %判断分数fs是否处在合格级别上 jb=' C '; %定义为C级 else %分数fs不处于以上任何级别上 jb=' D '; %定义为D级 end %选择结构结束 ► jb=dengji(81) %调用函数dengji ◄ jb =B 例4-3：Fibonacci 数组的元素满足Fibonacci规则： ： ， , . 求出该数组中第一个大于10000的元素。 n=100; %给定一个较大的n作为数列的位置 a=[1,1]; %设定数列的初始值 for i=3:n %从第3个元素开始循环递推生成后续元素 p=a(i-1)+a(i-2); %前两个元素之和生成后续元素p a=[a,p]; %将刚产生的元素p放置到数组a的最后，拼接成新的数组 if p>10000 %判断将刚产生的元素p是否超过10000 break; %跳出所在的for循环 end %选择结构结束 end %循环结构结束 disp([a]) %显示所生成的数列，最后一个元素a(length(a))为所求的元素 例4-4：动态显示数列极限 的逼近过程。 clear;clf; %清除内存变量，清理图形窗口 hold on %开启图形保持功能以便重复画点 axis([0,150,2,2.8]); %设置坐标窗口 grid %画出坐标网格 for n=2:2:150 %让n从2开始，建立循环 an=(1+1/n)^n; %计算数列的值 plot(n,an,'r.','markersize',15); %画出相应的坐标点，点的大小为15 pause(0.1); %暂停0.1秒后开始下一循环 fprintf('n=%d an=%.4f\n',n,an); %显示出每次结算结果（在命令窗口中） end %循环结构结束 问题1：对于数列 ， >0为常数，可以证明该数列收敛，且 。显然，这个结论提供了一个求平方根 的近似方法，试编制一个函数程序，对任意给定的正实数A，求出 的近似值（精确到 ）。 function [a1,k]=jixian(A) %定义函数jixian，输入常数A，返回极限近似值和迭代次数 if A<=0 %如果A不为正，则显示错误信息且a1=-1,k=-1并终止程序 disp('data error!'); a1=-1;k=-1; return %终止程序返回调用处 end a0=1; k=0; %给定数列的初始值a0和累计循环的变量k r=1; %设定存放相邻两次迭代值差的绝对值，赋初值1 while r>=1.0e-5 %当误差没达到 时继续循环 k=k+1; %循环次数累加一次 a1=(a0+A/a0)/2; %计算下一个元素a1 r=abs(a1-a0); %计算前后两次迭代值差的绝对值，存放在r中 a0=a1; %将a1赋给a0作为新的初值，迭代下一个元素 fprintf('k=%.0f a1=%.8f\n',k,a1) %显示每次迭代结果 end %循环结构结束 ► [a,b]=jixian(2) %调用函数jixian k=1 a1=1.50000000 k=2 a1=1.41666667 k=3 a1=1.41421569 k=4 a1=1.41421356 ◄ a = 1.4142 ◄ b = 4 问题2：对于任意一个正整数，都可以判断其是质数还是合数，这一点在一些有关数论问题中是经常用到的。但当一个正的奇数比较大时，手工来判断是否为质数往往不很容易。现在要求编制一个函数程序，对任意一个正整数，判断出它是质数还是合数，若是质数，则返回值1；若是合数，返回值0，同时给出两个因数；若输入非正数，则返回值-1，并提示错误。 function k=hezhishu(M) %定义函数hezhishu，输入正整数M，返回合质数的标志k k=1; %赋数k初值为1，默认为质数 if M<=1 disp(' data error!');k=-1; %输入的数非合质数时，输出错误信息 return %程序终止，返回 end if M>=4 %输入的数M从4开始，2，3为质数 for m=2: fix(sqrt(M)) %从2到 开始循环 if mod(M,m)==0 %M被m整除，即M为合数 k=0; %赋数k为0，表示M为合数 fprintf('%.0f=%.0f×%.0f\n',M,m,M/m) %显示M=m×n break %循环终止 end end end ► k=hezhishu(1259) %调用函数hezhishu ◄ k = 1 %表明1259为质数 ► k=hezhishu(1159) %调用函数hezhishu ◄ 1159 =19×61 ◄ k = 0 %表明1159为合数 ► k=hezhishu(-8) %调用函数hezhishu ◄ data error! ◄ k = -1 %表明输入数据错误 问题3：设某一建筑公司要筹建一批A、B、C三种类型的楼房，已知每栋楼房的投资和售价分别为：A类投资90万，售价115万；B类投资110万，售价150万；C类投资170万，售价205万。现在该公司有资金1250万，要求每类楼房至少建一栋，最多不超过5栋，那么如何设计建楼 方案 气瓶 现场处置方案 .pdf气瓶 现场处置方案 .doc见习基地管理方案.doc关于群访事件的化解方案建筑工地扬尘治理专项方案下载 ，在资金充分利用的前提下能获得最大利润？ clear;clc; t=[90,110,170]; %楼房的投资单价向量 p=[115,150,205]; %楼房的售价向量 z=1250; %总资本 r=p-t; %楼房的利润向量 D=[];k=0; %D为存放方案及利润的结果矩阵，k为行标 for a=1:5 % a为A类楼房的数量 for b=1:5 % b为B类楼房的数量 for c=1:5 % c为C类楼房的数量 f=[a,b,c]; % f为建楼方案即各类楼房数量向量 if f * t'<=z % 计算f方案投资额度，如果投资允许，方案有效 zr=f * r'; % 计算f方案的利润 k=k+1; %结果矩阵的行数加1 D(k,:)=[f , zr]; %方案及利润的结果存放矩阵D的k行 end end end end [R, I]=max(D( : ,4)); %求出矩阵D第4列（利润）中最大值所在的行I disp([D(I, : )]) %显示出矩阵D的第I行，前三列即为建楼最优方案 问题4：设 (0,0)为一导弹发射点，发现位于B(0,100)处一架敌机沿水平方向逃离（如图），随即发射一枚导弹予以打击，现已知导弹时刻对准敌机，且速率为飞机速率的两倍（设飞机速度为1）。试编程模拟导弹打击敌机的动态过程，并实时给出飞机和导弹的位置坐标。如果敌机飞行60单位距离之外即逃出我方空域，那么，要想在我方空域内击落敌机，则导弹的速度至少应提高到敌机速度的多少倍？ clear;clc;clf; hold on axis([0 100 0 120]); grid A=[0,0]; %标记导弹的初始位置 B=[0,100]; %标记敌机的初始位置 d=norm(A-B); %计算两者之间的距离 k=0; % k为循环统计量，即为迭代的次数 v=1;dt=1;K=2; %给定飞机速度v、时间间隔dt、速度关系系数K的值 while k<10000 %开始循环，次数不超过10000 k=k+1; %循环次数累加1 B=B+[v*dt,0]; %计算敌机位置 plot(A(1),A(2),'r.' ,'markersize',15); %标记导弹位置 plot(B(1),B(2),'b.' ,'markersize',15); %标记敌机位置 e=B-A;d=norm(e); %计算两者之间的距离 fprintf('k=%.0f B(%.0f,100) A(%.2f,%.2f) d=%.2f\n',k,B(1),A(1),A(2),d) if d<=0.5 %如果飞机与导弹距离很小，可认为击中目标 break %终止循环 end e=e/d; %计算导弹飞行方向向量 A=A+K*v*dt*e; %计算导弹位置 pause(0.2) %暂停0.2秒，放慢画图演示过程 end 实验5 开普勒方程近似解与方程求根 例 5-1 用“二分法”求方程 的近似根（误差 ）. f=inline(' x-0.5*sin(x)-1'); a=1; b=2; dlt=1.0e-5; k=1; while abs(b-a)>dlt c=(a+b)/2; if f(c)= =0 break; elseif f(c)*f(b)<0 a=c; else b=c; end fprintf('k=%d, x=%.5f\n',k,c); k=k+1; end 例 5-2 用“切线法”求方程 的近似根（误差 ）. f=inline(' x-0.5*sin(x)-1'); df=inline('1-0.5*cos(x)'); d2f=inline('0.5*sin(x)'); a=1; b=2; dlt=1.0e-5; if f(a)*d2f(a)>0 x0=a; else x0=b; end m=min(abs(df(a)),abs(df(b))); k=0; while abs(f(x0))>m*dlt k=k+1; x1=x0-f(x0)/df(x0); x0=x1; fprintf('k=%d x=%.5f\n',k,x0); end 例5-3 求方程组 的近似解. function f=group1(x) f=[sin(x(1))+x(2)+x(3)^2*exp(x(1))-4; x(1)+x(2)*x(3); x(1)*x(2)*x(3)+2]; ► [x,fval]=fsolve('group1',[1,1,1]) x = 1.4142 -1.3701 1.0322 fval = 1.0e-012 * 0.1155 0.0007 -0.0071 例5-4 求方程组 的近似解. function f=group2(x) f=[9*x(2)^2-12*x(1)-54*x(2)+61; x(1)*x(2)-2*x(1)+1]; ► [x,fval]=fsolve('group2',[0,0]) x = 1.0902 1.0828 fval = 1.0e-011 * 0.1044 -0.0324 ► [x,fval]=fsolve('group2',[-2,2]) x = -1.5682 2.6377 fval = 1.0e-007 * 0.3218 0.0260 实验6 Logistic方程求解与混沌 3 数值实验与分析 clc;clf; x=0.1; y=[ ]; r=1.2; %改变取值得到相应的图形 hold on axis([0 100 0 1]) for i=1:100 x=r*x*(1-x); y=[y,x]; plot(i,x,'k.','markersize',10) fprintf('x(%d)=%.10f\n',i,x); end t=1:100; plot(t,y,'k-'); grid 步4： clear;clf; axis([0,4,0,1]); grid; hold on for r=0:0.3:3.9 x=[0.1]; for i=2:150 x(i)=r*x(i-1)*(1-x(i-1)); end pause（0.5） for i=101:150 plot(r,x(i),'k.'); end text(r-0.1,max(x(101:150))+0.05,['\it{r}=',num2str(r)]) end 步6： clear;clf; hold on axis([2.7,4,0,1]); grid for r=2.7:0.005:3.9 x=[0.1]; for i=2:150 x(i)=r*x(i-1)*(1-x(i-1)); end pause(0.1) for i=101:150 plot(r,x(i),'k.'); end end 实验7 的计算与数值积分 2．级数逼近 clc;clear; n=0; r=1; p=0; k=-1; while r>=1.0e-5 n=n+1; k=k*(-1); p1=p+k/(2*n-1); r=abs(4*(p1-p)); fprintf('n=%.0f,p=%.10f\n',n,4*p1); p=p1; end MATLAB程序： clear; n=0; r=1; p=0; k=-1; a=1; b=1; while r>=1.0e-5 n=n+1; k=k*(-1); a=4*a;b=9*b; p1=p+k/(2*n-1)*(2/a+3/b); r=abs(4*(p1-p)); fprintf('n=%.0f,p=%.10f\n',n,4*p1); p=p1; end 3．数值积分 （1）梯形法 f=inline('4./(1+x.*x)'); x=0:0.1:1; y=f(x); p=trapz(x,y); fprintf('p=%.6f\n',p) 再对分点加密，精度有所提高： f=inline('4./(1+x.*x)'); x=0:0.01:1; y=f(x); p=trapz(x,y); fprintf('p=%.6f\n',p) （2）复化梯形法 clear; f=inline('4./(1+x.*x)'); a=0;b=1; n=1; h=(b-a)/n; t1=h/2*(f(a)+f(b)); er=1;k=1; while er>1.0e-5 s=0; for i=1:n s=s+f(a+(i-1/2)*h); end t2=(t1+h*s)/2; er=abs(t2-t1); fprintf('n=%.0f,p=%.6f,r=%.6f\n',k,t2,er); n=2*n; h=h/2; t1=t2; k=k+1; end （3）抛物线（Simpson）方法 clear; f=inline('4./(1+x.*x)'); a=0;b=1; n=1; z=quad(f,a,b); fprintf('z=%.10f\n',z) 实验8 河流流量估计与数据插值 （1）多项式插值 function p=lagrange(x,y) %输入数据点坐标向量x,y，输出插值多项式系数p L=length(x); A=ones(L); for j=2:L A(:,j)=A(:,j-1).*x'; end X=inv(A)*y'; for i=1:L p(i)=X(L-i+1); end 例8-1 已知观测数据 x 1 2 3 4 5 y -1 1.5 2.1 3.6 4.9 求其插值多项式曲线。 x=[1 2 3 4 5]; y=[-1 1.5 2.1 3.6 4.9]; plot(x,y,'k.','markersize',15) axis([1 5 -1 5]) grid; hold on p=lagrange(x,y); t=1:0.1:5; u=polyval(p,t); plot(t,u,'r-') 例8-2 已知观测数据 x 0 .1 .2 .3 .4 .5 .6 .7 .8 .9 1 y -.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2 求其插值多项式曲线。 x=0:0.1:1; y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]; plot(x,y,'k.','markersize',15) axis([0 1 -2 16]) grid; hold on p=lagrange(x,y); t=0:0.01:1; u=polyval(p,t); plot(t,u,'r-') 例8-3 对函数 ，在[-5,5]上以1为步长进行划分作Lagrange插值，观察函数曲线（虚线）与插值曲线（实线）的变化。 x=-5:0.1:5; y=1./(1+20*x .*x); plot(x,y,'k--','linewidth',2) axis([-5 5 -1.2 6]) grid; hold on x=-5:5; y=1./(1+20*x .*x); p=lagrange(x,y); t=-5:0.1:5; f=polyval(p,t); plot(t,f,'r-') （2）分段线性插值 x=0:0.1:1; y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]; plot(x,y,'k.','markersize',15) axis([0 1 -2 16]) grid; hold on t=0:0.01:1; u=interp1(x,y,t); plot(t,u,'r-') （3） 三次样条插值 x=0:0.1:1; y=[-.447 1.978 3.28 6.16 7.08 7.34 7.66 9.56 9.48 9.3 11.2]; plot(x,y,'k.','markersize',15) axis([0 1 -2 16]) grid;pause（0.5） hold on t=0:0.01:1; u=spline(x,y,t); plot(t,u,'r-') （1）画出河床观测点的散点图 clf;clear x=0:5:100; y=[0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 ... 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0]; y1=10-y; plot(x,y1,'k.','markersize',18); axis([0 100 0 10]); grid （2）利用分段线性插值绘制河床曲线 clf;clear x=0:5:100; y=[0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 ... 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0]; y1=10-y; plot(x,y1,'k.','markersize',15); axis([0 100 2 10]) grid;hold on t=0:100; u=interp1(x,y1,t); plot(t,u) S=100*10-trapz(x,y1); p=sqrt(diff(x).^2+diff(y1).^2); L=sum(p); fprintf('S=%.2f , L=%.2f\n',S,L) （3）利用样条插值绘制河床曲线 x=0:5:100; y=[0 2.41 2.96 2.15 2.65 3.12 4.23 5.12 6.21 5.68 4.22 ... 3.91 3.26 2.85 2.35 3.02 3.63 4.12 3.46 2.08 0]; y1=10-y; plot(x,y1,'k.','markersize',15); axis([0 100 2 10]) grid;hold on t=0:100; u=spline(x,y1,t); plot(t,u) S=100*10-trapz(t,u); p=sqrt(diff(t).^2+diff(u).^2); L=sum(p); fprintf('S=%.2f , L=%.2f\n',S,L) 实验9 人口预测与数据拟合 例9-1 已知观测数据 x 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 y -.447 1.