MATLAB软件与基础数学实验.doc

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}表示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 采用不同形式（直角坐标、参数、极坐标），画出单位圆的图形。
（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