软件与基础数学实验

1、MATLAB 软件与基础数学实验 实验 1 MATLAB 基本特性与基本运算 例 1-1 求12+2(7-4)32的算术运算结果。 （12+2*（7-4） ）/32 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 设 75,24ba ，计算 |)tan(| |)|sin(| ba ba 的值。 a=pi/180*(-

2、24); %转换为弧度值且不显示 b=pi/180*75; %转换为弧度值且不显示 z=sin(abs(a)+abs(b)/sqrt(tan(abs(a+b) %计算结果并显示 z =0.8888 例 1-6 设三角形三边长为 2, 3, 4cba ，求此三角形的面积。 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 设 101 654 321 A ， 112 311 021 B ，计算 | ,AABBA ， 1 A。 A=1,

3、2,3;4,5,6;1,0,1; B=-1 2 0;1 1 3;2 1 1; BA ans= 0 4 3 5 6 9 3 1 2 BA* Ans=7 7 9 13 19 21 1 3 1 )(det A % det 为求方阵的行列式命令 ans = -6 )(inv A %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，只要输入下列语句：

4、 A(2,3)= -1； 例 1-9 分别画出函数 xxycos 2 和x x z sin 在区间-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) % 画函数

5、y 的曲线,默认为蓝色（参看实验 2） grid %在第一个子窗中加坐标网格 subplot(1,2,2) %在第二个子窗中绘图 plot(x,z,linewidth,2) % 画函数 z 的曲线,默认为蓝色（参看实验 2） grid %在第二个子窗中加坐标网格 例 1-10 试求方程组 4 3 2 201 624 121 X 的解。 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.00

6、00 1.5000 5.0000 1.0000 -1.0000 -3.0000 x=c*b %矩阵左逆乘，结果为方程组的解 x = -30.0000 22.5000 -13.0000 X=ab %用除法直接求方程组的解 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 试求矩阵方程 111

7、 321 201 624 121 X 的解。 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 建立同时计算 n bay)( 1 ， n bay)( 2 的函数。即任给 a,b,n 三个数，返 回 y1,y2. function y1 , y2=fun1(a , b , n) % fun1 is a function used by DEMO y1=(a+b)n,

8、 y2=(a-b)n % Copyright by XJTU y1=(a+b).n ; y2=(a-b).n; 例 1-13 设 22 11 ( )6 (0.3)0.01(0.9)0.04 f x xx ，试画出在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

9、)-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 的极限 ans=cos(x) 例 1-18：设 )sin(),(yyxyxf n ，求 ., 2 2 2 yx f y f y f x f syms x y n %

10、声明符号变量，注意变量间必须用空格分开 fx=xn*y+sin(y); %建立符号函数 diff(fx) %对变量 x(默认)求一阶导数（偏导数） ans =xn*n/x*y 即 ynxn 1 diff(fx, y) %对变量 y 求一阶导数（偏导数） ans =xn+cos(y) diff(fx, y, 2) %对变量 y 求二阶导数（偏导数） ans =-sin(y) diff(diff(fx,x), y) %先对 x 求导再对 y 求导（二阶混合偏导数） ans = xn*n/x 即 1n nx 例 1-19：求 dx x xy 2 1， dy x xy t 0 2 1， dy x xy

11、 dx x 0 2 1 0 1， .)( 1 0 1 0 1 0 yxx dzzyxdydx syms x y z %声明符号变量，注意变量间必须用空格分开 f1=x*y/(1+x2) ; %建立符号函数 f2=x+y+z; int(f1) %对 f1 关于变量 x(默认)求不定积分 ans =1/2*y*log(1+x2) %即 )1ln( 2 1 2 xy syms t int(f1,0, t ) %对 f1 关于变量 x(默认)在0,t上求定积分 ans =1/2*log(1+t2)*y %即 )1ln( 2 1 2 ty int(int(f1,y,0, sqrt(x),x,0,1 )

12、%对 f1 先求对 y 的积分再求对 x 的积分（二重积 分） ans =1/2-1/8*pi %即 . 8 1 2 1 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) %求级数 k 1 3 1 2 1 1 (ans=inf 即) symsum(1/(k*(k+1),1,inf) %求级数 ) 1( 1 32 1 21 1 kk (ans=1) symsum(a*1/3k,k,0,inf) %求

13、级数 k aaa a 333 2 (ans= 3/2*a) 泰勒展开（taylor） syms x fy=1/(1+x+x2) f=taylor(fy) %求 fx 对自变量 x(默认)在 x=0 点(默认)泰勒展开前 6 项(默认) f=taylor(fy,8,1) %求 fx 对自变量 x(默认)在 x=1 点泰勒展开式前 8 项 方程求根（solve） fx=sym(a*x2+b*x+c) ; %建立符号函数 solve(fx) %求方程 fx=0 的符号解 ans = 1/2/a*(-b+(b2-4*a*c)(1/2) 1/2/a*(-b-(b2-4*a*c)(1/2) syms b

14、solve(fx, b ) %求方程 fx=0 关于变量 b 的符号解 ans = -(a*x2+c)/x 微分方程(组)求解（dsolve） dsolve(Dy=5) %求方程 y=5 的通解，默认自变量为 t ans = 5*t+C1 dsolve(Dy=x, x) %求方程 y=x 的通解，指定自变量为 x ans =1/2*x2+C1 dsolve(D2y=1+Dy, y(0)=1, Dy(0)=0) %求方程 y=1+y满足 y(0)=1,y(0)=0 的特解 ans = -t+exp(t) 即 t ety x,y=dsolve(Dx=x+y,Dy=2*x) %求方程组 2xy yx

15、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( 3/1)2( 3/1 )2( 3/1)( 3/1 2 2121 2 2121 tt tt ecceccy ecceccx 实验 2 MATLAB 绘制二维、三维图形 例 2-1 在子图形窗口中画出 2 , 0 上正弦、余弦曲线。 x=0:0.1*pi:2*pi; %按步长赋值生成 x 向量 y=

16、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 画出 2 , 0 上正弦、余弦曲线并对线型加粗、点型加大，重新定置坐标系以 及加注相关说明和注释

17、。 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 itx) %标记横坐标轴, itx表示 x 为斜体 ylabel(Variable ity) %标记纵坐标轴 title(Sine and Cosine Cruves) %标记图名 text(2.5,

18、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 分别在两个图形窗口画出填充一正方形和极坐标方程 2cos2sin2r 的图形。 h1=figure; %打开第一个图形窗口，返回其图标识号（句柄）h1 x=0 1 1 0 0; %闭合图形的顶点横坐标向量 y=0 0 1 1 0; %闭

19、合图形的顶点纵坐标向量 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,

20、linewidth,3) %对第二个窗口中曲线加粗 例 2-4 在-2.5,2.5上画出函数 2 x ey 的直方图和阶梯图。 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

21、 Plot of a Bell Curve ) %添加图形标题 例 2-5 采用不同形式（直角坐标、参数、极坐标） ，画出单位圆 1 22 yx 的图形。 （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); %生成

22、单位圆上的函数值向量 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, 10, 0t 上一段曲线。 t=0:pi/50:10*pi; %生成参数 t 数组 X=sin(t); %生成螺旋线 X 数组 Y=cos(t); %生成螺旋线 Y 数组 Z=t; %生成螺旋线 Z 数组 plot3(X,Y,Z

23、, k-, linewidth,3) %画螺旋线 grid 例 2-7 画出矩形域-1,1 -1,1上旋转抛物面： 22 yxz 。 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; %生成 22 yxz 函数值矩阵 subplot(1,2,1) mesh(X,Y,Z) ; %在第一个子图中画 22 yxz 网格曲面 subplot(1,2,2) surf(X,Y,Z) ; %在第二个子图中画 22 yxz 光滑曲面 shading f

24、lat ; %对曲面 22 yxz 平滑并除去网格 例 2-8 在圆形域 1 22 yx 上绘制旋转抛物面： 22 yxz 。 x=linspace(-1,1,300); %分割-1,1区间生成 x y=x; %生成 y X,Y=meshgrid(x,y); %生成矩形域-1,1X-1,1网格节点坐标矩阵 Z=X.2+Y.2; %生成 22 yxz 函数值矩阵 i=find(Z1); %找出圆域 1 22 yx 之外的函数值（z1）坐标点 i Z(i)=NaN; %对圆域 1 22 yx 之外的坐标点 i 处函数值进行 “赋空” subplot(1,2,1) mesh(X,Y,Z) ; %在第

25、一个子图中画 22 yxz 网格曲面 subplot(1,2,2) surf(X,Y,Z) ; %在第二个子图中画 22 yxz 光滑曲面 shading flat ; %对曲面 22 yxz 平滑并除去网格 例 2-9 画出 22 22 sin yx yx z 在 5 . 7| , 5 . 7|yx 上的图形。 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 有一组实验数据如下表所示，试绘图表示。 时

26、间1 2 3 4 5 6 7 8 9 数据 112.51 13.54 15.60 15.92 20.64 24.53 30.24 50.00 36.34 数据 29.87 20.54 32.21 40.50 48.31 64.51 72.32 85.98 89.77 数据 310.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

27、 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 xlabel(time);ylabel(data); axis(0 10 0 100); text(6.5,25.5,leftarrowdata1); % leftarrow 表示画一左箭头,且在标识前 text(3,43.8,data2rightarrow); % rightarrow 表示画一右箭头,且在标识后

28、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=%.0fn,i,s) %逐行显示出累加求和的过程 end %循环结构结束 例

29、 3-2：求正整数 n 的阶乘：p=12 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=%.0fn,i,p) %逐行显示出 i! end %循环结构结束 例 3-3：根据麦克劳林公式可以得到 e1+1+1/2!+1/3!+1/n!，试求 e 的近似值。 clear;clc; %清除内存变

30、量，清理命令窗口 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=%.8fn,i,s) %逐行显示出第 i 次 e 的近似值 end %循环结构结束 例 3-4：对于数列 , 2 , 1,nn ，求其前 n 项和不超过 1000 时的 n 的值及和. cle

31、ar;clc; %清除内存变量，清理命令窗口 n=0; %设定正整数并赋初值 0 s=0; %设定存放累加和的变量 s 并赋初值 0 while s=1.0e-8 %当近似值的精度 r 没达到 9 10时继续循环 k=k+1; %累计循环次数并作为下一个正整数 k p=p*k; %计算 k 的阶乘 p r=1/p; %计算前后两次近似值的误差 r s=s+r; %计算 e 的近似值 s fprintf(k=%.0f, s=%.10fn,k,s) %逐行显示出第 k 次 e 的近似值 s end %循环结构结束 实验 4 MATLAB 选择结构与应用实验 例 4-1：求任意有限数组 a=a(1)

32、,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; %

33、任意给定一数组 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 %分数

34、 fs 不处于以上任何级别上 jb= D ; %定义为 D 级 end %选择结构结束 jb=dengji(81) %调用函数 dengji jb =B 例 4-3：Fibonacci 数组的元素满足 Fibonacci 规则： n a ： 1 21 aa ， 12 kkk aaa , , 3 , 2 , 1k . 求出该数组中第一个大于 10000 的元素。 n=100; %给定一个较大的 n 作为数列的位置 a=1,1; %设定数列的初始值 for i=3:n %从第 3 个元素开始循环递推生成后续元素 p=a(i-1)+a(i-2); %前两个元素之和生成后续元素 p a=a,p; %将

35、刚产生的元素 p 放置到数组 a 的最后，拼接成新的数组 if p10000 %判断将刚产生的元素 p 是否超过 10000 break; %跳出所在的 for 循环 end %选择结构结束 end %循环结构结束 disp(a) %显示所生成的数列，最后一个元素 a(length(a)为所求的元素 例 4-4：动态显示数列极限 )( ) 1 1 (ne n a n n 的逼近过程。 clear;clf; %清除内存变量，清理图形窗口 hold on %开启图形保持功能以便重复画点 axis(0,150,2,2.8); %设置坐标窗口 grid %画出坐标网格 for n=2:2:150 %让

36、 n 从 2 开始，建立循环 an=(1+1/n)n; %计算数列的值 plot(n,an,r.,markersize,15); %画出相应的坐标点，点的大小为 15 pause(0.1); %暂停 0.1 秒后开始下一循环 fprintf(n=%d an=%.4fn,n,an); %显示出每次结算结果（在命令窗口中） end %循环结构结束 问题 1：对于数列 )( 2 1 : 1 n nnn a A aaa ), 2 , 1 , 0(n ， Aa, 0 0 0 为常数，可 以证明该数列收敛，且 Aan n lim 。显然，这个结论提供了一个求平方根 A的近似方 法，试编制一个函数程序，对任

37、意给定的正实数 A，求出 A的近似值（精确到 5 10） 。 function a1,k=jixian(A) %定义函数 jixian，输入常数 A，返回极限近似值和迭代 次数 if A=1.0e-5 %当误差没达到 6 10时继续循环 k=k+1; %循环次数累加一次 a1=(a0+A/a0)/2; %计算下一个元素 a1 r=abs(a1-a0); %计算前后两次迭代值差的绝对值，存放在 r 中 a0=a1; %将 a1 赋给 a0 作为新的初值，迭代下一个元素 fprintf(k=%.0f a1=%.8fn,k,a1) %显示每次迭代结果 end %循环结构结束 a,b=jixian(2

38、) %调用函数 jixian k=1 a1=1. k=2 a1=1. k=3 a1=1. k=4 a1=1. a = 1.4142 b = 4 问题 2：对于任意一个正整数，都可以判断其是质数还是合数，这一点在一些有关数 论问题中是经常用到的。但当一个正的奇数比较大时，手工来判断是否为质数往往不很容 易。现在要求编制一个函数程序，对任意一个正整数，判断出它是质数还是合数，若是质 数，则返回值 1；若是合数，返回值 0，同时给出两个因数；若输入非正数，则返回值-1， 并提示错误。 function k=hezhishu(M) %定义函数 hezhishu，输入正整数 M，返回合质数的标 志 k

39、k=1; %赋数 k 初值为 1，默认为质数 if M=4 %输入的数 M 从 4 开始，2，3 为质数 for m=2: fix(sqrt(M) %从 2 到 M 开始循环 if mod(M,m)=0 %M 被 m 整除，即 M 为合数 k=0; %赋数 k 为 0，表示 M 为合数 fprintf(%.0f=%.0f%.0fn,M,m,M/m) %显示 M=mn break %循环终止 end end end k=hezhishu(1259) %调用函数 hezhishu k = 1 %表明 1259 为质数 k=hezhishu(1159) %调用函数 hezhishu 1159 =19

40、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 栋，那么如何设计建楼方案，在资金充分利用的前提下能获得最大利润？ clear;clc; t=90,110,170; %楼房的投资单价向量 p=1

41、15,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 行

42、 end end end end R, I=max(D( : ,4); %求出矩阵 D 第 4 列（利润）中最大值所在的行 I disp(D(I, : ) %显示出矩阵 D 的第 I 行，前三列即为建楼最优方 案 问题 4：设 0 A (0,0)为一导弹发射点，发现位于 B(0,100)处一架敌机沿水平方向逃离 （如图） ，随即发射一枚导弹予以打击，现已知导弹时刻对准敌机，且速率为飞机速率的两 倍（设飞机速度为 1） 。试编程模拟导弹打击敌机的动态过程，并实时给出飞机和导弹的位 置坐标。如果敌机飞行 60 单位距离之外即逃出我方空域，那么，要想在我方空域内击落敌 机，则导弹的速度至少应提高到敌

43、机速度的多少倍？ 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 km*dlt k=k+1; x1=x0-f(x0)/df(x0); x0=x1; fprintf(k=%d x=%.5fn,k,x0); end 例 5-3 求方程组 1 2 123 123 123 sin()40

44、0 2 x xxx e xx x x x x 的近似解. 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 求方程组 2 212 121 91254610 210 xxx x xx 的近似解. function f=group2(x) f=9*x(2)2-12*x(1)-

45、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*

46、x*(1-x); y=y,x; plot(i,x,k.,markersize,10) fprintf(x(%d)=%.10fn,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,itr=,num2str

47、(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=%.10fn,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=%.10fn,n,4*p1); p=p1; end 3数值积分 （1）梯形法 f=inline(4./(1+x.*x); x=0:0.1:1; y=f(x); p