数学模型与实验

1、数学模型与实验（一）Matlab简介函数的数值表达初等函数的数值直接表示Sqrt sin cos log sin(2) ans = 0.9093 log(3) ans = 1.0986向量的生成1.直接输入向量，格式要求是：向量元素需要用“ ” 括起来，元素之间可以用空格、逗号或分号分隔； a=1,2,3 a = 1 2 3 sin(a)ans = 0.8415 0.9093 0.14112.利用冒号表达式生成向量冒号表达式的基本形式为x=x0：step：xnx0表示向量的首元素值，xn表示向量的尾元素数值限，step表示从第二个元素开始，元素数值大小与前一个元素之大小的差值。 a=1:2:1

2、2 a = 1 3 5 7 9 11向量的基本运算1. 加（减）与数加（减） a-1 %这里的a即上面生成的a ans = 0 2 4 6 8 102. 数乘 a*2 ans = 2 6 10 14 18 223. 点积 dot（a，b）返回向量a和b的数量点积。a和b必须同维。 a=1 2 3 b=3 4 5 dot(a,b)ans = 264.叉积cross（a，b）返回向量a和b的叉积向量。a，b必须为三维向量。 a=1 2 3 b=3 4 5 cross(a,b)ans = -2 4 -2矩阵的生成1.直接输入小矩阵输入矩阵时要以” “为标识，同行元素之间可用空格或”，“分隔，行与行之

3、间要用”；“或回车符分隔，若不想看到中间结果，可以用”；“结束。 a=1 2 3;1 ,1 1;4,5 6 a = 1 2 3 1 1 1 4 5 6 2.直接导入文件矩阵的运算矩阵的加和减 a=1 2 3;1 ,1 1;4,5 6; b=3 4 5;1 2 1;2 3 1; c=a+bc = 4 6 8 2 3 2 6 8 7矩阵的乘法 e=b,5 5 5e = 3 4 5 5 1 2 1 5 2 3 1 5 f=a*ef = 11 17 10 30 6 9 7 15 29 44 31 75矩阵的逆运算【例】求下面A矩阵的逆。解： a= 2 1 -3 -1；3 1 0 7；-1 2 4 -2

4、； 1 0 -1 5；解： a=2 1 -3 -1；3 1 0 7；-1 2 4 -2； 1 0 -1 5； inv(a) ans = -0.0471 0.5882 -0.2706 -0.9412 0.3882 -0.3529 0.4824 0.7647 -0.2235 0.2941 -0.0353 -0.4706 -0.0353 -0.0588 0.0471 0.2941 2.矩阵的行列式运算 a1=det(a)a1 = -853.秩函数 rank(a)ans = 4特征值函数 x,y=eig(a)% x为特征向量矩阵，y为特征值矩阵x = -0.5843 0.9223 -0.1387 -

5、0.2449i -0.1387 + 0.2449i 0.7160 0.0904 -0.7828 -0.7828 -0.3806 -0.0066 -0.0803 + 0.3491i -0.0803 - 0.3491i 0.0333 -0.3756 -0.4212 - 0.0489i -0.4212 + 0.0489i y = -1.1228 0 0 0 0 2.5266 0 0 0 0 5.2981 + 1.3755i 0 0 0 0 5.2981 - 1.3755i几种常用的工具阵zeros（M,N）生成M*N阶的全0阵zeros(size(A)生成与A同阶的全0阵eye（M,N）生成M*N阶

6、的单位阵ones（M,N）生成M*N阶的全1阵rand（M,N）生成M*N阶的随机矩阵数组的运算数组的加减法与矩阵的加减法运算完全相同，符号统一“+”，“-” a= 2 1 -3 -1;3 1 0 7;-1 2 4 -2;1 0 -1 5; a+1 ans =3 2 -2 0 4 2 1 8 0 3 5 -1 2 1 0 6数组的乘除法是指两同维数组间对应元素之间的乘除法，它的运算符为“.*”和“./”或“.” a=1 2 3;2 3 4;3 4 5; b=1 1 1;2 2 2;3 3 3; a./bans = 1.0000 2.0000 3.0000 1.0000 1.5000 2.000

7、0 1.0000 1.3333 1.6667数组的幂运算符为“.”,他表示每个数组元素的幂运算，同矩阵的幂运算符不同。 a3 a.3ans = ans = 132 192 252 1 8 27 192 279 366 8 27 64 252 366 480 27 64 125多项式的表示方法对于多项式用以下的行向量表示这样就把多项式问题转化为向量问题【例如】输入多项式 p=1 -5 6 -33; poly2sym(p)ans = x3 - 5*x2 + 6*x 33%其中poly2sym是符号工具箱中的函数，可将多项式向量表示为符号多项式的形式多项式的拟合polyfit (X, Y，n) 其中

8、X、Y为拟合数据 ，n为拟合多项式的阶数【例】用5阶多项式对0,/2上的正弦函数进行最小二乘拟合。 x=0:pi/20:pi/2; y=sin(x); a=polyfit(x,y,5); x1=0:pi/30:pi*2; y1=sin(x1);y2=a(1)*x1.5+a(2)*x1.4+a(3)*x1.3+a(4)*x1.2+a(5)*x1+a(6); plot(x1,y1,b-,x1,y2,r*) legend(原曲线,拟合曲线) axis(0,7,-1.2,4)函数拟合观察数据（x，y）；对于线性拟合模型y=ax+b：1.利用A=polyfit(x,y,1)得到系数向量 其中A(1)=a

9、，A(2)=b 所得线性模型为 y=A(1)*x+A(2)2.解线性方程组 对于线性拟合模型y=ax+b，等价于 令A=(a，b), 则 那么MatLab图形绘制功能基本绘图函数命 令含 义plot建立向量或矩阵各队队向量的图形loglogx、y轴都取对数标度建立图形semilogxx轴用于对数标度，y轴线性标度绘制图形semilogyy轴用于对数标度，x轴线性标度绘制图形title给图形加标题xlabel给x轴加标记ylabel给y轴加标记text在图形指定的位置上加文本字符串gtext在鼠标的位置上加文本字符串grid打开网格线x=0:0.001:10; % 0到10的1000个点的x座标

10、 y=sin(x); % 对应的y座标 plot(x,y); % 绘图（可以设定显示范围）画椭圆a = 0:pi/50:2*pi; X = cos(a)*3; %参数方程Y = sin(a)*2;plot(X,Y);xlabel(x), ylabel(y);title(椭圆) 线条和颜色控制符:.点y黄。小圆圈m棕色xx标记c青色+加号r红色-实线g绿色*星号b蓝色：虚线w白色-.点划线h黑色条形图当资料点数量不多时，长条图是很适合的表示方式： close all; % 关闭所有的图形视窗 x=1:10; y=rand(size(x); bar(x,y) 误差图x = linspace(0,2

11、*pi,30); y = sin(x); e = std(y)*ones(size(x); % std:标准差errorbar(x,y,e) 极坐标图形 theta=linspace(0, 2*pi); r=cos(4*theta); polar(theta, r); x=linspace(0,10,50); y=sin(x).*exp(-x/3); stairs(x,y); %阶梯图x=linspace(0,10,50); y=sin(x).*exp(-x/3); stem(x,y); %针状图符号表达式画图在许多的场合，将表达式可视化是有利的。MATLAB提供了函数ezplot来完成该任务

12、。 y= 16*x2+64*x+96 y=16*x2+64*x+96 ezplot(y)三维曲线和曲面 1空间曲线绘制函数plot3( )格式一：plot3(x,y,z)功能：若x、y、z为相同长度的向量，则根据向量x、y、z绘制空间三维曲线。如果x、y、z为同阶矩阵，则绘制对应列的多条曲线。格式二：plot3(x, y, z, s)功能：按字符串s设置的线型、颜色、标记符号绘制三维空间曲线。字符串s的设置与plot命令的设置相同。 2三维网线图函数mesh、meshc和meshz MATAL提供的函数mesh用来绘制三维网线图。 格式一：mesh(X,Y,Z) 功能：根据矩阵X、Y和Z绘制彩

13、色的空间三维网线图。X、Y和Z中对应的元素为三维 空间上的点，点与点之间用线连接。其中网线的颜色随着网点高度的改变而改变。 格式二：mesh(x,y,Z) 功能：n维向量x、m维向量y和mn矩阵Z绘制网线图，节点的坐标为 (x(j),y(i),Z(i,j)，网线的颜色随着网点高度的改变而改变。三维曲面图函数surf及surfc格式一：surf(X, Y, Z, C) 功能：绘制由四个矩阵所指定的带色参数的网状表面图。视角是由view所指定。轴的刻度决定于X，Y及Z的范围，或当前对轴的设定；颜色范围由C指定。格式二：surf(X，Y，Z)功能：将C设为与Z相等，则颜色与网的高度成正比。三维螺旋线

14、 t=0:pi/50:10*pi; plot3(sin(t),cos(t),t) grid 空间中的曲线t=linspace(0,20*pi, 501); plot3(t.*sin(t), t.*cos(t), t); %注意用点乘 .*t=linspace(0, 10*pi, 501); plot3(t.*sin(t), t.*cos(t), t, t.*sin(t), t.*cos(t), -t); 立体网状图 x=linspace(-2, 2, 25); % 在x轴上取25点 y=linspace(-2, 2, 25); % 在y轴上取25点 xx,yy=meshgrid(x, y);

15、% xx和yy都是21x21的矩阵 zz=xx.*exp(-xx.2-yy.2); % 计算函数值，zz也是21x21的矩阵 mesh(xx, yy, zz); % 画出立体网状图 网状图填补x=linspace(-2, 2, 25); % 在x轴上取25点 y=linspace(-2, 2, 25); % 在y轴上取25点 xx,yy=meshgrid(x, y); % xx和yy都是21x21的矩阵 zz=xx.*exp(-xx.2-yy.2); % 计算函数值，zz也是21x21的 矩阵 surf(xx, yy, zz); % 画出立体曲面图 peaks函数peaks函数，可产生一个凹凸

16、有致的曲面，包含了三个局部极大点及三个局部极小点，其方程式为： x,y,z=peaks; waterfall(x,y,z); 等高线图meshc同时画出网状图与等高线： x,y,z=peaks; meshc(x,y,z); surfc同时画出曲面图与等高线：x,y,z=peaks; surfc(x,y,z); 剔透玲珑球X0,Y0,Z0=sphere(30); %产生单位球面的三维坐标X=2*X0;Y=2*Y0;Z=2*Z0; %产生半径为2的球面的三维坐标surf(X0,Y0,Z0); %画单位球面shading interp %采用插补明暗处理hold on; mesh(X,Y,Z);ho

17、ld off %画外球面hidden off %产生透视效果axis off %不显示坐标轴彗星状轨迹图 【*例】简单二维示例。（请指令窗中运行以下指令）n=10;t=n*pi*(0:0.0005:1);x=sin(t);y=cos(t);plot(x,y,g); hold on;comet(x,y,0.02);hold off卫星返回地球的运动轨线示意 R0=1; %以地球半径为一个单位a=12*R0;b=9*R0;T0=2*pi;%T0是轨道周期T=5*T0;dt=pi/100;t=0:dt:T;f=sqrt(a2-b2);%地球与另一焦点的距离th=12.5*pi/180;%卫星轨道与x

18、-y平面的倾角E=exp(-t/20);%轨道收缩率x=E.*(a*cos(t)-f);y=E.*(b*cos(th)*sin(t);z=E.*(b*sin(th)*sin(t);plot3(x,y,z,g)%画全程轨线X,Y,Z=sphere(30);X=R0*X;Y=R0*Y;Z=R0*Z; %获得单位球坐标grid on,hold on,surf(X,Y,Z),shading interp%画地球x1=-18*R0;x2=6*R0;y1=-12*R0;y2=12*R0;z1=-6*R0;z2=6*R0;axis(x1 x2 y1 y2 z1 z2) %确定坐标范围view(117 37)

19、,comet3(x,y,z,0.02),hold off%设视角、画运动轨线 求解线性方程组 MATLAB语句实现： a=2,3,-1；8,2,3；45,3,9； b=2；4；23； x=inv(a)*b求多项式方程的根roots(p) solve(p)例. 求方程 的所有根 p=1 -4 9 -10 r=roots(p) 或 s1=sym(x3-4*x2+9*x-10); solve(s1)求方程组的根和超越方程解x,y= solve(x+y=1,x-11*y=5,x,y)求超越方程 的根，其中p，r为常数ff=sym(p*sin(x)=r) solve(ff,x)求极限syms x; y1

20、=(1+4*x)(1/x); y2=(exp(x)-1) /x; limit(y1,x,0) limit(y2,x,0)求导数与偏导数 syms x; f=3*x3+5*x+1; diff(f,1) syms x y; %对x偏导数 z=x2*sin(2*y); B=diff(z,x)求解不定积分syms x b;f=cos(b*x);f1=int(f,x)f1 = 1/b*sin(b*x)ezplot(f1)求定积分多重积分 syms x ; y=(x2+sin(x)/(1+x2); int(y,x,-1,1) syms x y ; int(int(x2+y2+1,y,x,x+1),x,0,

21、1)求微分方程通解和特解 dsolve(Dy=1+y2, x) dsolve(Dy=1+y2, y(0)=1 ,x)dsolve(x2*D2y+x*Dy+(x2-1/2)*y=0,y(pi/2)=2,Dy(pi/2)=-2/pi,x)MATLAB语言编程介绍MATLAB中各种命令可以完成许多单一的任务，对于某些较为复杂的问题，仅靠现有的命令或函数来解决，往往是难以达到目的 。为此，要运用MATLAB编程语言编制程序，形成M-文件。程序是使计算机完成各项运算的命令集，运行一个编制好的程序，计算机会从第一条命令行开始，一行接一行地执行相应的命令，直到终止。程序一般包含：数据输入，数据处理和结果输出

22、三部分，其中数据处理是核心。程序编写调试完成后，需要存盘，形成永久性文件，可以随时对它进行调用或修改。文件 名以字母开头，但不能用专用变量名，如pi等。 MATLAB中每一个命令都是一个M-文件。MATLAB程序类型脚本文件函数文件function y=mean(x)% MEAN Average or mean Value% For vectors,Mean(x) is the mean value of X% For matrices, Mean(x) is a row vector containing% the mean value of each columnm,n=size(x)if

23、 m=1 m=n;endy=sum(x)/m;函数定义行H1行函数帮助文件函数体M-文件中输入、输出命令 直接赋值输入： 提示对话输入(input命令) x=2; x=input(请输入参数 x=); a=1,2;3,4; a=input(请输入矩阵 a=); s= any string!; s=input(Please input s=); 直接输出： 格式控制输出(fprintf命令) x fprintf(x=%.0f, y=%.5fn,pi,pi); disp(a,b); fprintf(x=%5g, y=%10.5gn,N,pi);M-文件中循环控制命令(for命令)格式： for i

24、=n1:(step):n2 commands; end 作用：重复执行命令集commands.例2：求奇数和：s=1+3+5+(2k-1) n=input(请输入 n=); s=0; for i=1:2:n s=s+i; fprintf(i=%.0f, s=%.0fn,i,s); end例1：求和：s=1+2+n n=input(请输入 n=); s=0; for i=1:n s=s+i; fprintf(i=%.0f, s=%.0fn,i,s); end例3：求阶乘：p=12 3 n=n! n=input(请输入 n= ); p=1; for i=1:n p=p*i; fprintf( i=

25、%.0f, p=%.0fn ,i,p); end例4：求e：e=1+1+1/2!+1/3!+1/n! n=input(请输入 n= ); p=1;e=1; for i=1:n p=p*i; p1=1/p; e=e+p1; fprintf( i=%.0f, p=%.0f, e=%.8f n ,i,p,e); end例5：求pi：pi/4=1-1/3+1/5-1/7+1/9-1/11+ n=input(请输入 n= ); pi=0;k=-1; for i=1:2:n k=(-1)*k; pi=pi+k*1/i; fprintf( i=%.0f, pi=%.8f n ,i,4*pi); endM-文

26、件中条件循环命令(while命令)格式： while (condition is true) commands; end 作用：重复执行命令集commands.例6：求和：s=1+2+n3000 s=0;k=1; while s3000 s=s+k; fprintf( k=%.0f, s=%.0fn ,k,s); k=k+1; end例7：求e：e=1+1+1/2!+1/3!+ (1/n!=1.0e-8 p=p*i; r=1/p; e=e+r; fprintf( i=%.0f, r=%.9f, e=%.9f n ,i,r,e); i=i+1; endM-文件中选择控制命令(if命令)单项选择控

27、制格式： if (condition is true) commands; end 作用：若条件成立，则执行命令集 commands. 否则，不执行。例8：求n个实数中最大的数M. a=input(请输入数组 an= ); m,n=size(a); M=a(1); for i=2:n if M0 sn=1;elseif x=0 sn=0;else sn=-1;endfprintf(x=%.5f, sn=%.0fn,x,sn);Matlab编程 (续)程序设计技术 循环向量化计算模式 目的: 向量计算替代循环模式 实例: i初始化 tic tic for t=0:0.01:100 t=0.0:0.01:100 i=i+1; y(i)=sin(t) y=sin(t) end toc toc 耗时7.58 耗时0.01 Matlab编程 (续)程序设计技术(续) 阵列预分配空间 目的: 避免循环过程空间分配的耗时 实例: y(k)=0.75y(k-1)-0.125y(k-2)+2u(k) clear all,tic T=0.001; t=0:T:16 u=sin(2*pi*t); % y=zeros(size(u) y(1)=2*u(1); y(2)=0.75*y(1)+2*u(2);耗时0.562(%)/0.0160for k=3:fix(16