《科学计算与matlab语言》,中国大学MOOC, 中南大学
https://www.icourse163.org/learn/CSU-1002475002#/learn/content
使用中的技巧发现 Pwd 输出当前工作目录,支持Linux命令。
diary - 将命令行窗口文本记录到日志文件中
save - 将工作区变量保存到文件中
load - 将文件变量加载到工作区中
0 课程导入 求x^2-3x+1=0的根。
方法一:利用MATLAB多项式求根函数roots来求根。
绘图:
1 2 3 4 x =-5 :0 .1 :5 ;y1 =x.*x-3 *x+1 ;y2 =zeros(size(x));plot (x, y1 , x, y2 );
方法二 : 利用求单变量非线性方程根的函数fzero,求方程在某个初始点附近的实根。
1 2 3 4 f =@(x) x*x-3 *x+1 x1 =fzero(f, 0.5 ) x2 =fzero(f, 2.5 )
方法三:利用最优化工具箱中的方程求根函数fsolve。
1 2 3 4 f =@(x) x*x-3 *x+1 x1 =fsolve(f, 0.5 , optimset('Display' , 'off' )) x2 =fsolve(f, 2.5 , optimset('Display' , 'off' ))
方法四:利用solve函数求方程的符号解,即求得的解是一个表达式。
1 2 3 4 syms x x = solve(x ^2 -3 *x +1 ) x = eval(x )
1 matlab基础知识 1.1 matlab系统环境 1.2 matlab数值数据 例1 分别求一个三位正整数的个位数字、十位数字和百位数字。
1 2 3 4 m=345 ; m1=rem (m,10 ) m2=rem (fix (m/10 ),10 ) m3=fix (m/100 )
例2 求[1,100]区间的所有素数。
1 2 3 4 x=1 :100 ; k=isprime (x); k1=find (k); p=x(k1)
1.3 变量及其操作 例1 计算表达式的值,并将结果赋给变量z,然后显示计算结果。
1 2 3 x=sqrt (7 )-2 i ; y=exp (pi /2 ); z=(5 +cosd (47 ))/(1 +abs (x-y))
1.4 MATLAB矩阵的表示方法 1.5 矩阵元素的引用 1.6 MATLAB基本运算 例1 当x=0.1、0.4、0.7、1时,分别求y=sin x cos x的值。
1 2 x=0.1 :0.3 :1 ; y=sin (x).*cos (x)
例2 建立3阶方阵A,判断A的元素是否为偶数。
1 2 A =[24 ,35 ,13 ;22 ,63 ,23 ;39 ,47 ,80 ] P=rem (A,2 )==0
例3 水仙花数是指各位数字的立方之和等于该数本身的三位正整数。求全部水仙花数。
1 2 3 4 5 6 m=100 :999 ; m1=rem (m,10 ); m2=rem (fix (m/10 ),10 ); m3=fix (m/100 ); k=find (m==m1.*m1.*m1+m2.*m2.*m2+m3.*m3.*m3) s=m(k)
1.7 字符串处理 例1 建立一个字符串向量,然后对该向量做如下处理: ① 取第1~5个字符组成的子字符串。 ② 将字符串倒过来重新排列。 ③ 将字符串中的小写字母变成相应的大写字母,其余字符不变。 ④ 统计字符串中小写字母的个数。
1 2 3 4 5 6 ch='ABc123d4e56Fg9' ; subch=ch(1 :5 ) revch=ch(end :-1 :1 ) k=find (ch>='a' &ch<='z' ) ch(k)=ch(k)-('a' -'A' ) length (k)
专题一总结 2 matlab 矩阵处理 2.1 特殊矩阵 例1 产生5阶两位随机整数矩阵A,再产生均值为0.6、方差为0.1的5阶正态分布随机矩阵B,验证(A+B)I=IA+BI(I为单位矩阵)。
1 2 3 4 A=fix (10 +(99 -10 +1 )*rand (5 )); B=0.6 +sqrt (0.1 )*randn (5 ); C=eye (5 ); (A+B)*C==C*A+B*C
例2 产生8阶魔方阵,求其每行每列元素的和。
1 2 3 M=magic (8 ); sum(M(1 ,:)) sum(M(:,1 ))
例3 生成5阶帕斯卡矩阵,验证它的逆矩阵的所有元素也为整数。
1 2 3 format rat P=pascal (5 ) inv(P)
2.2 矩阵变换 例1 先建立5×5矩阵A,然后将A的第一行元素乘以1,第二行乘以2,…,第五行乘以5。
1 2 3 A=[7 ,0 ,1 ,0 ,5 ;3 ,5 ,7 ,4 ,1 ;4 ,0 ,3 ,0 ,2 ;1 ,1 ,9 ,2 ,3 ;1 ,8 ,5 ,2 ,9 ] D=diag (1 :5 ); D*A
要将A的各列元素分别乘以对角阵的对角线元素。
1 2 3 A=[7 ,0 ,1 ,0 ,5 ;3 ,5 ,7 ,4 ,1 ;4 ,0 ,3 ,0 ,2 ;1 ,1 ,9 ,2 ,3 ;1 ,8 ,5 ,2 ,9 ] D=diag (1 :5 ); A*D
例2 验证魔方阵的主对角线、副对角线元素之和相等。
1 2 3 4 5 6 7 A=magic (5 ); D1=diag (A); sum(D1) B=flipud (A); D2=diag (B); sum(D2)
例3 用求逆矩阵的方法解线性方程组Ax=b。
1 2 3 4 A=[1 ,2 ,3 ;1 ,4 ,9 ;1 ,8 ,27 ]; b=[5 ;-2 ;6 ]; x=inv(A)*b x=A\b
2.3 矩阵求值 例1 验证det(A-1)=1/det(A)。
1 2 3 4 format rat A=[1 ,3 ,2 ;-3 ,2 ,1 ;4 ,1 ,2 ] det(inv(A)) 1 /det(A)
例2 求3~20阶魔方阵的秩。
1 2 3 4 5 6 7 for n=3 :20 r(n)=rank(magic (n)); end bar(r) grid on axis([2 ,21 ,0 ,20 ]) [3 :20 ;r(3 :20 )]
例3 求2~10阶希尔伯特矩阵的条件数。
1 2 3 4 5 for n=2 :10 c(n)=cond(hilb (n)); end format long c'
2.4 矩阵的特征值与特征向量 例1
1 2 3 4 5 6 R=[-1 ,2 ,0 ;2 ,-4 ,1 ;1 ,1 ,-6 ]; S=[1 ,2 ;2 ,3 ]; A=[R,zeros (3 ,2 );zeros (2 ,3 ),S]; [X1,d1]=eig(R) [X2,d2]=eig(S) [X3,d3]=eig(A)
例2
1 2 3 4 5 6 7 x=[0 ,0.5 ,0.5 ,3 ,5.5 ,5.5 ,6 ,6 ,3 ,0 ;0 ,0 ,6 ,0 ,6 ,0 ,0 ,8 ,1 ,8 ]; A=[1 ,0.5 ;0 ,1 ]; y=A*x; subplot(2 ,2 ,1 ); fill(x(1 ,:),x(2 ,:),'r' ); subplot(2 ,2 ,2 ); fill(y(1 ,:),y(2 ,:),'r' );
2.5 稀疏矩阵 求下列三对角线性方程组的解。
1 2 3 4 5 6 7 8 kf1=[1 ;1 ;2 ;1 ;0 ]; k0=[2 ;4 ;6 ;6 ;1 ]; k1=[0 ;3 ;1 ;4 ;2 ]; B=[kf1,k0,k1]; d=[-1 ;0 ;1 ]; A=spdiags(B,d,5 ,5 ); f=[0 ;3 ;2 ;1 ;5 ]; x=A\f
专题二总结 3 matlab 程序流程控制 3.1 顺序结构程序 例1 有一线段AB,A的坐标为(1,1),B的坐标为(4.5,4.5),求AB的长度,以及黄金分割点C的坐标。
1 2 3 4 5 6 a=input('a=' ); b=input('b=' ); c=a+0.618 *(b-a); s=abs (a-b); disp (s)disp (c)
3.2 用if语句实现选择结构 例1 输入一个整数,若为奇数则输出其平方根,否则输出其立方根。
1 2 3 4 5 6 7 x=input('请输入x的值:' ); if rem (x,2 )==1 y=sqrt (x); else y=x^(1 /3 ); end y
例2 输入一个字符,若为大写字母,则输出其对应的小写字母;若为小写字母,则输出其对应的大写字母;若为数字字符则输出其对应数的平方,若为其他字符则原样输出。
1 2 3 4 5 6 7 8 9 10 c=input('请输入一个字符:' ,'s' ); if c>='A' && c<='Z' disp (lower(c)) elseif c>='a' && c<='z' disp (upper(c)) elseif c>='0' && c<='9' disp (str2double(c)^2 ) else disp (c) end
3.3 用switch语句实现选择结构 例1 输入一个英文单词,判断它是否以元音字母开头。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 c=input('请输入一个单词:' ,'s' ); switch c(1 ) case {'A' ,'E' ,'I' ,'O' ,'U' ,'a' ,'e' ,'i' ,'o' ,'u' } disp ([c,'以元音字母开头' ]); otherwise disp ([c,'以辅音字母开头' ]); end c=input('请输入一个单词:' ,'s' ); if findstr(c(1 ),'AEIOUaeiou' )>0 disp ([c,'以元音字母开头' ]); else disp ([c,'以辅音字母开头' ]); end
例2 PM2.5是指大气中直径小于或等于2.5微米的可入肺颗粒物,是衡量空气质量的重要指标。假定空气质量等级以PM2.5数值划分为6级。 PM2.5数值[0,35)空气质量为优,[35,75)为良,[75,115)为轻度污染,[115,150)为中度污染,[150,250)为重度污染,大于等于250为严重污染。编写程序,输入PM2.5数值,输出空气质量等级。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 g=input('请输入PM2.5值:' ); switch fix (g) case num2cell (0 :34 ) disp ('空气质量优' ); case num2cell (35 :74 ) disp ('空气质量良好' ); case num2cell (75 :114 ) disp ('空气质量轻度污染' ); case num2cell (115 :149 ) disp ('空气质量中度污染' ); case num2cell (150 :249 ) disp ('空气质量重度污染' ); otherwise disp ('空气质量严重污染' ); end
3.4 用for语句实现循环 计算圆周率π。
(1)利用无穷级数展开式求π的近似值。
1 2 3 4 5 6 7 8 y=0 ; g=-1 ; n=input('n=?' ); for i =1 :n g=-g; y=y+g/(2 *i -1 ); end pai=4 *y
用向量求和的方法实现程序:
1 2 3 4 n=input('n=?' ); x=1 :2 :(2 *n-1 ); y=(-1 ).^(2 :n+1 )./x; pai=sum(y)*4
(2)利用定积分的近似值求π的近似值。
1 2 3 4 5 6 7 8 9 10 11 12 a=0 ; b=1 ; n=input('n=?' ); h=(b-a)/n; x=a:h:b; f=sqrt (1 -x.*x); s=[]; for k=1 :n s1=(f(k)+f(k+1 ))*h/2 ; s=[s,s1]; end pai=4 *sum(s)
(3)利用蒙特卡洛法求π的近似值。
1 2 3 4 5 6 7 8 9 10 11 12 > s=0 ; > n=input('n=?' ); > for i =1 :n > x=rand (1 ); > y=rand (1 ); > if x*x+y*y<=1 > s=s+1 ; > end > end > pai=s/n*4
3.5 用while语句实现循环 例1 从键盘输入若干个数,当输入0时结束输入,求这些数的平均值和它们之和。
1 2 3 4 5 6 7 8 9 10 11 12 msum=0 ; n=0 ; x=input('Enter a number (end in 0):' ); while x~=0 msum=msum+x; n=n+1 ; x=input('Enter a number (end in 0):' ); end if n>0 msum mean =msum/n end
例2 求[100,200]之间第一个能被21整除的整数。
1 2 3 4 5 6 7 for n=100 :200 if rem (n,21 )~=0 continue end n break end
例3 用筛选法求某自然数范围内的全部素数。
1 2 3 4 5 6 7 8 9 10 m=input('m=' ); p=1 :m; p(1 )=0 ; for i =2 :sqrt (m) for j =2 *i :i :m p(j )=0 ; end end n=find (p~=0 ); p(n)
3.6 函数文件的定义与调用 例1 编写函数文件,求半径为r的圆的面积和周长。
1 2 3 function [s,p] =fcircle (r) s=pi *r*r; p=2 *pi *r;
例2 已知=(40)/((30)+(20)) ①当()=+10 ln(^2+5)时,的值是多少。 ②当()=1×2+2×3+3×4+⋯+×(+1)时,的值是多少。 分别用匿名函数和函数文件定义函数()。 第②问的函数文件f2.m。
1 2 3 4 5 function f =f2 (n) f=0 ; for k=1 :n f=f+k*(k+1 ); end
脚本文件mf.m。
1 2 3 f1=@(n) n+10 *log (n*n+5 ); y1=f1(40 )/(f1(30 )+f1(20 )) y2=f2(40 )/(f2(30 )+f2(20 ))
3.7 函数的递归调用 例1 利用函数的递归调用,求n!。
函数文件fact.m如下:
1 2 3 4 5 6 function f =fact (n) if n<=1 f=1 ; else f=fact (n-1 )*n; end
在脚本文件a.m中调用函数文件fact.m,求n!。
1 2 3 n=input('Please input n=' ); s=fact (n); disp (s)
例2 Fibonacci数列定义如下: f1=1 f2=1 fn=fn-1+fn-2 (n>2) 编写递归调用函数求Fibonacci数列的第n项,然后调用该函数验证Fibonacci数列的如下性质:
首先建立函数文件ffib.m。
1 2 3 4 5 6 function f =ffib (n) if n>2 f=ffib(n-1 )+ffib(n-2 ); else f=1 ; end
建立程序文件test.m。
1 2 3 4 5 6 F=[]; for k=1 :20 F=[F,ffib(k)*ffib(k)]; end sum(F) ffib(20 )*ffib(21 )
3.8 函数参数与变量作用范围 例1 利用nargin和nargout建立函数文件test.m。
1 2 3 4 5 6 7 8 function fout =test (a,b,c) if nargin==1 fout=a; elseif nargin==2 fout=a+b; elseif nargin==3 fout=(a*b*c)/2 ; end
例2 利用全局变量建立函数文件wad.m。
1 2 3 function f =wad (x,y) global ALPHA BETAf=ALPHA*x+BETA*y;
专题三总结 4 matlab 绘图 4.1 二维曲线 例1 绘制一条折线。
1 2 3 x=[2.5 , 3.5 , 4 , 5 ]; y=[1.5 , 2.0 , 1 , 1.5 ]; plot (x, y)
例2 绘制sin(x)、sin(2x)、sin(x/2)的函数曲线。
1 2 3 4 5 6 7 8 9 x=linspace (0 ,2 *pi ,100 ); y=[sin (x); sin (2 *x); sin (0.5 *x)]; plot (x,y) t=0 :0.01 :2 *pi ; t1=t'; x=[t1, t1, t1]; y=[sin (t1), sin (2 *t1), sin (0.5 *t1)]; plot (x,y)
例3 采用不同个数的数据点绘制正弦函数曲线,观察曲线形态。
1 2 3 4 t1=linspace (0 , 2 *pi , 10 ); t2=linspace (0 , 2 *pi , 20 ); t3=linspace (0 , 2 *pi , 100 ); plot (t1, sin (t1), t2,sin (t2)+1 , t3, sin (t3)+2 )
例4 用不同线型和颜色在同一坐标内绘制曲线及其包络线。
1 2 3 4 5 6 x=(0 :pi /50 :2 *pi )'; y1=2 *exp (-0.5 *x)*[1 ,-1 ]; y2=2 *exp (-0.5 *x).*sin (2 *pi *x); x1=0 :0.5 :6 ; y3=2 *exp (-0.5 *x1).*sin (2 *pi *x1); plot (x, y1, 'k:' , x, y2, 'b--' , x1, y3, 'rp' )
例5 绘制函数sin(1/)的图形。
1 2 3 x=0 :0.005 :0.2 ; y=sin (1. /x); plot (x,y)
例6 采用fplot函数绘制函数sin(1/x)。
1 fplot(@(x) sin (1. /x),[0 ,0.2 ], 'b' )
例7 已知螺旋线的参数方程,绘制曲线。
1 fplot(@(t)t.*sin (t), @(t)t.*cos (t), [0 ,10 *pi ], 'r' )
4.2 绘制图形的辅助标注 例1 绘制[-2π,2π]区间的正弦曲线并给图形添加标题。
1 2 3 4 5 x=-2 *pi :0.05 :2 *pi ; y=sin (x); plot (x,y)title('y=sin(x)' ) title({'MATLAB' , 'y=sin(x)' })
\bf 加粗
\it 斜体
\rm 正体
Color(’k’),fontsize(11)
绘制[-2π,2π]区间的正弦曲线并给x轴添加标签。
1 2 3 4 5 x=-2 *pi :0.05 :2 *pi ; y=sin (x); plot (x,y)title('y=sin(x)' ) xlabel('-2\pi \leq x \leq 2\pi' )
在前面的图形中添加文字说明。Text(x,y,txt);gtext(txt)
1 2 text(-2 *pi , 0 , '-2{\pi}' ) text(3 , 0.28 , '\leftarrow sin(x)' )
例2 绘制不同频率的正弦曲线并用图例标注曲线。Legend()
1 2 3 x = linspace (0 , 2 *pi , 100 ); plot (x, [sin (x); sin (2 *x); sin (3 *x)])legend ('sin(x)' , 'sin(2x)' , 'sin(3x)' )
绘制边长为1的正方形,并设置坐标轴。Axis()
1 2 3 4 5 x = [0 , 1 , 1 , 0 , 0 ]; y = [0 , 0 , 1 , 1 , 0 ]; plot (x,y) axis([-0.1 , 1.1 , -0.1 , 1.1 ]) axis equal;
axis equal
axis square
axis auto
axis off
axis on
grid on
grid off
grid
box on
box off
box
例3 绘制sin(x)、sin(2x)、sin(x/2)的函数曲线并添加图形标注。
1 2 3 4 5 6 7 8 9 10 11 x=linspace (0 ,2 *pi ,100 ); y=[sin (x); sin (2 *x); sin (0.5 *x)]; plot (x,y)axis([0 7 -1.2 , 1.2 ]) title('不同频率正弦函数曲线' ); xlabel('Variable X' ); ylabel('Variable Y' ); text(2.5 , sin (2.5 ), 'sin(x)' ); text(1.5 , sin (2 *1.5 ), 'sin(2x)' ); text(5.5 , sin (0.5 *5.5 ), 'sin(0.5x)' ); legend ('sin(x)' ,'sin(2x)' ,'sin(0.5x)' )grid on
例4 用图形保持功能绘制两个同心圆。Hold on ;hold off;hold
1 2 3 4 5 6 7 8 t = linspace (0 ,2 *pi ,100 ); x = sin (t); y = cos (t); plot (x, y, 'b' )hold on; plot (2 *x, 2 *y, 'r--' )grid on axis([-2.2 2.2 -2.2 2.2 ]) axis equal
划分2×2子图 subplot(m,n,p)
1 2 3 4 5 6 subplot(2 ,2 ,1 ); x=linspace (0 ,2 *pi ,60 ); y=sin (x); plot (x,y);title('sin(x)' ); axis ([0 ,2 *pi ,-1 ,1 ]);
划分多子图
1 2 3 4 5 6 7 8 9 10 11 12 13 x=linspace (0 ,2 *pi ,60 ); subplot(2 ,2 ,1 ) plot (x,sin (x)-1 );title('sin(x)-1' );axis ([0 ,2 *pi ,-2 ,0 ]) subplot(2 ,1 ,2 ) plot (x,cos (x)+1 );title('cos(x)+1' );axis ([0 ,2 *pi ,0 ,2 ]) subplot(4 ,4 ,3 ) plot (x,tan (x));title('tan(x)' );axis ([0 ,2 *pi ,-40 ,40 ]) subplot(4 ,4 ,8 ) plot (x,cot (x));title('cot(x)' );axis ([0 ,2 *pi ,-35 ,35 ])
4.3 其他形式的二维图形 例1 绘制1/的直角线性坐标图和三种对数坐标图。
对数坐标图semilogx();seilogy();loglog()
1 2 3 4 5 6 7 8 9 10 11 12 13 14 x=0 :0.1 :10 ; y=1. /x; subplot(2 ,2 ,1 ) plot (x,y) title('plot(x,y)' );grid on subplot(2 ,2 ,2 ) semilogx(x,y) title('semilogx(x,y)' );grid on subplot(2 ,2 ,3 ) semilogy(x,y) title('semilogy(x,y)' );grid on subplot(2 ,2 ,4 ) loglog(x,y) title('loglog(x,y)' );grid on
例2 按极坐标方程ρ=1-sin t绘制心形曲线。Polar(theta,rho,options)
1 2 3 4 5 6 7 8 t = 0 :pi /100 :2 *pi ; r = 1 -sin (t); subplot(1 ,2 ,1 ) polar(t,r) subplot(1 ,2 ,2 ) t1 = t-pi /2 ; r1 = 1 -sin (t1); polar(t,r1)
例3 绘制分组条形图。bar(y,style);barh();stacked:堆积分组;groud:簇状分组(默认)
1 2 3 4 5 6 7 y=[1 ,2 ,3 ,4 ,5 ; 1 ,2 ,1 ,2 ,1 ; 5 ,4 ,3 ,2 ,1 ]; subplot(1 ,2 ,1 ) bar(y) title('Group' ) subplot(1 ,2 ,2 ) bar(y, 'stacked' ) title('Stack' )
例4 下表是某公司2015~2017年家电类商品1月份的销售数据,绘制条形图对比数据。
1 2 3 4 5 6 x=[2015 ,2016 ,2017 ]; y=[68 ,80 ,115 ,98 ,102 ; 75 ,88 ,102 ,99 ,110 ;81 ,86 ,125 ,105 ,115 ];bar(x, y) title('Group' );
例5 绘制服从高斯分布的直方图。hist();直角坐标系;rose():极坐标系。
1 2 3 4 5 6 7 8 y=randn (500 ,1 ); subplot(2 ,1 ,1 ); hist(y); title('高斯分布直方图' ); subplot(2 ,1 ,2 ); x=-3 :0.2 :3 ; hist(y,x); title('指定区间中心点的直方图' )')
例6 绘制高斯分布数据在极坐标下的直方图。
1 2 3 4 y=randn (500 ,1 ); theta=y*pi ; rose(theta) title('在极坐标下的直方图' )
例7 某次考试优秀、良好、中等、及格、不及格的人数分别为:7、17、23、9、4,试用扇形统计图作成绩统计分析。Pie函数:扇形图;area函数:面积图
1 2 3 4 5 score = [5 , 17 , 23 , 9 , 4 ]; ex = [0 ,0 ,0 ,0 ,1 ]; pie(score, ex) legend ('优秀' , '良好' , '中等’, ' 及格', ' 不及格', … ' location', 'eastoutside' )
例8 以散点图形式绘制桃心曲线。Scatter散点图;stairs阶梯图;stem杆图
1 2 3 4 t = 0 :pi /50 :2 *pi ; x = 16 *sin (t).^3 ; y = 13 *cos (t)-5 *cos (2 *t)-2 *cos (3 *t)-cos (4 *t); scatter (x,y,'rd' ,'filled' )
例9 已知向量A、B,求A+B,并用矢量图表示。Compass 罗盘图;feather羽毛图;quiver箭头图
1 2 3 4 5 6 7 8 A=[4 ,5 ]; B=[-10 ,0 ]; C=A+B; hold on;quiver(0 , 0 , A(1 ), A(2 )); quiver(0 , 0 , B(1 ), B(2 )); quiver(0 , 0 , C(1 ), C(2 )); text(A(1 ),A(2 ),'A' );text(B(1 ),B(2 ),'B' ); text(C(1 ),C(2 ),'C' ); axis ([-12 , 6 , -1 , 6 ]) grid on
4.4 三维曲线 例1 绘制一条空间折线。Plot3(x,y,z);fplot3
1 2 3 4 5 6 x=[0.2 , 1.8 , 2.5 ]; y=[1.3 , 2.8 , 1.1 ]; z=[0.4 , 1.2 , 1.6 ]; plot3 (x, y, z)grid on axis([0 , 3 , 1 , 3 , 0 , 2 ]);
例2 绘制螺旋线
1 2 3 4 5 6 7 8 9 10 t=linspace (0 , 10 *pi , 200 ); x=sin (t)+t.*cos (t); y=cos (t)-t.*sin (t); z=t; subplot(1 , 2 , 1 ) plot3 (x, y, z)grid on subplot(1 , 2 , 2 ) plot3 (x(1 :4 :200 ), y(1 :4 :200 ), z(1 :4 :200 ))grid on
例3 在空间不同位置绘制5条正弦曲线。
1 2 3 4 5 6 7 8 9 10 11 12 t=0 :0.01 :2 *pi ; t=t'; x=[t, t, t, t, t]; y=[sin (t), sin (t)+1 , sin (t)+2 , sin (t)+3 , sin (t)+4 ]; z=x; plot3 (x,y,z)这个例子也可以采用以下代码实现。 t=0 :0.01 :2 *pi ; x=t; y=[sin (t); sin (t)+1 ; sin (t)+2 ; sin (t)+3 ; sin (t)+4 ]; z=x; plot3 (x,y,z)
例4 绘制三条不同长度的正弦曲线。
1 2 3 4 5 t1=0 :0.01 :1.5 *pi ; t2=0 :0.01 :2 *pi ; t3=0 :0.01 :3 *pi ; plot3 (t1,sin (t1),t1, t2,sin (t2)+1 ,t2, …t3,sin (t3)+2 ,t3)
例5 绘制空间曲线
1 2 3 4 5 6 7 t=0 :pi /50 :6 *pi ; x=cos (t); y=sin (t); z=2 *t; plot3 (x,y,z,'p' )xlabel('X' ),ylabel('Y' ),zlabel('Z' ); grid on
例6 绘制墨西哥帽顶曲线。参数方程fplot3(funx,funy,funz,tlims[-5,5])
1 2 3 4 xt = @(t) exp (-t/10 ).*sin (5 *t); yt = @(t) exp (-t/10 ).*cos (5 *t); zt = @(t) t; fplot3(xt, yt, zt, [-12 , 12 ])
用红色点划线绘制墨西哥帽顶曲线。
1 2 3 4 xt = @(t) exp (-t/10 ).*sin (5 *t); yt = @(t) exp (-t/10 ).*cos (5 *t); zt = @(t) t; fplot3(xt, yt, zt, [-12 , 12 ], 'r-.' )
4.5 三维曲面: mesh三维网格图;surf三维曲面图;fmesh;fsurf;;meshgrid;meshc;meshz;surfc;surfl;
例1 绘制空间曲线。
1 2 3 4 5 6 x = 2 :6 ; y = (3 :8 )'; [X, Y] = meshgrid (x, y); Z = randn (size (X)); plot3 (X,Y,Z)grid on;
例2 绘制三维曲面图。
1 2 3 4 5 6 7 8 9 10 t = -2 :0.2 :2 ; [X, Y] = meshgrid (t); Z = X .* exp (-X.^2 - Y.^2 ); subplot(1 ,3 ,1 ) mesh(X,Y,Z); subplot(1 ,3 ,2 ) surf(X,Y,Z); subplot(1 ,3 ,3 ) plot3 (X,Y,Z); grid on
例3 用4种方式绘制函数$z=(x-1)^2+(y-2)^2-1$曲面图。 其中,x∈[0,2],y∈[1,3]。
1 2 3 4 5 6 7 8 9 10 [x,y]=meshgrid (0 :0.1 :2 ,1 :0.1 :3 ); z=(x-1 ).^2 +(y-2 ).^2 -1 ; subplot(2 ,2 ,1 ); meshc(x,y,z);title('meshc(x,y,z)' ) subplot(2 ,2 ,2 ); meshz(x,y,z);title('meshz(x,y,z)' ) subplot(2 ,2 ,3 ); surfc(x,y,z);title('surfc(x,y,z)' ) subplot(2 ,2 ,4 ); surfl(x,y,z); title('surfl(x,y,z)' )
例4 用cylinder函数分别绘制柱面、花瓶和圆锥面。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Sphere(n=20 );cylinder(R,n=20 ) subplot(1 ,3 ,1 ); [x,y,z]=cylinder; surf(x,y,z); subplot(1 ,3 ,2 ); t=linspace (0 ,2 *pi ,40 ); [x,y,z]= cylinder(2 +cos (t),30 ); surf(x,y,z); subplot(1 ,3 ,3 ); [x,y,z]= cylinder(0 :0.2 :2 ,30 ); surf(x,y,z); Peaks()
例5 用cylinder函数绘制两个相互垂直且直径相等的圆柱面的相交图形。
1 2 3 4 5 6 [x ,y,z]= cylinder(1,60); z=[-1z(2,:)]; surf(x ,y,z) hold on surf(y,z,x ) axis equal
例6 绘制螺旋曲面。
1 2 3 4 5 6 7 funx = @(u,v) u.*sin(v);funy = @(u,v) -u.*cos(v);funz = @(u,v) v;fsurf (funx,funy,funz,[-5 5 -5 -2 ])hold on fmesh (funx,funy,funz,[-5 5 -2 2 ])hold off
4.6 图形修饰处理
视点处理:方位角;仰角;view(az=-37.5,el=30)
View(x,y,z)
色彩处理[RGB]
colormap
剪裁处理
例1 绘制函数$z=(x-1)^2+(y-2)^2-1$曲面,并从不同视点展示曲面。
1 2 3 4 5 6 7 8 9 10 [x,y]=meshgrid (0 :0.1 :2 , 1 :0.1 :3 ); z=(x-1 ).^2 +(y-2 ).^2 -1 ; subplot(2 ,2 ,1 ); mesh(x,y,z) title('方位角=-37.5{\circ},仰角=30{\circ}' ) subplot(2 ,2 ,2 ); mesh(x,y,z) view(0 ,90 );title('方位角=0{\circ},仰角=90{\circ}' ) subplot(2 ,2 ,3 ); mesh(x,y,z) view(90 ,0 ); title('方位角=90{\circ},仰角=0{\circ}' ) subplot(2 ,2 ,4 ); mesh(x,y,z) view(-45 ,-60 ); title('方位角=-45{\circ},仰角=-60{\circ}' )
例2 创建一个灰色系列色图矩阵。
1 2 3 4 5 6 7 8 c = [0 , 0.2 , 0.4 , 0.6 , 0.8 , 1.0 ]'; cmap = [c, c, c]; surf(peaks) colormap(cmap) 或 cmap=gray(6 ); surf(peaks) colormap(cmap)
例3 使用同一色图,以不同着色方式绘制圆锥体。Shading faceted;shading flat;shading interp
1 2 3 4 5 6 7 8 [x,y,z]= cylinder(pi :-pi /5 :0 ,10 ); colormap(lines); subplot(1 ,3 ,1 ); surf(x,y,z); shading flat subplot(1 ,3 ,2 ); surf(x,y,z); shading interp subplot(1 ,3 ,3 ); surf(x,y,z);
例4 绘制3/4圆。将图形中需要裁减部分对应的函数值设置成NaN。
1 2 3 4 5 6 7 8 9 t = linspace (0 ,2 *pi ,100 ); x = sin (t); y = cos (t); p = y > 0.5 ; y(p)= NaN; plot (x,y)axis([-1.1 ,1.1 ,-1.1 ,1.1 ]) axis square grid on
例5 绘制3/4球面。
1 2 3 4 5 6 7 [X, Y, Z] = sphere(60 ); p = Z>0.5 ; Z(p) = NaN; surf(X, Y, Z) axis([-1 , 1 , -1 , 1 , -1 , 1 ]) axis equal view(-45 , 20 )
4.7 交互式绘图工具
“绘图”选项卡
图形窗口绘图工具:plottools;按钮;;图形选项板(新子图面板、变量面板、注释面板)、绘图浏览器、属性编辑器
图形窗口菜单和工具栏(工具栏、图形窗口菜单)
生成例1所需变量。
1 2 3 4 5 6 7 8 x=linspace (0 ,2 *pi ,100 ); y=sin (x); y1=sin (x); y2=sin (0.5 *x); y3=sin (2 *x); [u,v]=meshgrid (0 :0.1 :2 ,1 :0.1 :3 ); h=(u-1 ).^2 +(v-2 ).^2 -1 ; [X,Y,Z]= cylinder(0 :0.2 :2 ,30 );
专题四总结 5、数据分析与多项式计算 5.1 数据统计分析(matlab以列优先) Max/min/ 最大值、最小值
Mean/median 均值、中值
Sum/prod 求和、求积
Cumsum/cumprod 累加和、累乘积
Std/corrcoef 标准差、相关系数
样本标准差、总体标准差
Sort 排序
例1 求向量x的最大元素,其中x=[-43,72,9,16,23,47]。
1 2 3 4 \>> x = [-43 , 72 , 9 , 16 , 23 , 47 ] \>> y= max (x ) \>> [y, k]= max (x )
例2 求矩阵A(见课件)的每行及每列的最大元素,并求整个矩阵的最大元素。
1 2 3 4 5 6 \>> A= [13 ,-56 ,78 ;25 ,63 ,-235 ;78 ,25 ,563 ;1 ,0 ,-1 ]; \>> max (A) \>> max (A,[],2 ) \>> max (max (A))
例3 某学生宿舍的5位同学月生活费如向量x所示,其中,小明同学家境一般,请问他应该按什么标准向父母主张生活费额度才较为合理。x=[1200,800,1500,1000,5000]
1 2 3 4 \>> x=[1200 ,800 ,1500 ,1000 ,5000 ]; \>> mean (x) \>> median(x)
例4 求向量X=[1,2,3,4,5,6,7,8,9,10]的积与累乘积。
1 2 3 4 \>> X=[1 ,2 ,3 ,4 ,5 ,6 ,7 ,8 ,9 ,10 ]; \>> y1=prod(X) \>> y2=cumprod(X)
例5 生成满足正态分布的50000*4随机矩阵,用不同的形式求其各列之间的标准差。
1 2 3 4 5 6 7 8 9 10 11 \>> x=randn (50000 ,4 ); \>> y1=std(x,0 ,1 ) \>> y2=std(x,1 ,1 ) \>> x1=x'; \>> y3=std(x1,0 ,2 ); \>> y3' \>> y4=std(x1,1 ,2 ); \>> y4'
例6 某新产品上市,在上市之前,公司物流部门把新产品分配到不同地区的10个仓库进行销售。产品上市一个月后,公司要对各种不同的分配方案进行评估,以便在下一次新产品上市时进行更准确的分配,避免由于分配不当而产生的积压和断货。下表(见课件)是相关数据,请判断那种分配方案最为合理。
1 2 3 4 5 6 \>> A=[5032 ,6000 ,5100 ,5200 ;6532 ,6500 ,6600 ,5800 ; 5500 ,7000 ,5400 ,4800 ;4530 ,4000 ,4300 ,4200 ; 2300 ,2000 ,2200 ,2500 ;3254 ,3000 ,3500 ,3000 ; 8095 ,9000 ,7800 ,8500 ;7530 ,8000 ,7000 ,7500 ; 3841 ,3200 ,3500 ,3200 ;4500 ,5200 ,4800 ,4000 ]; \>> corrcoef(A)
例7 对下列矩阵(见课件)做各种排序。
1 2 3 4 5 6 \>> A=[1 ,-8 ,5 ;4 ,12 ,6 ;13 ,7 ,-13 ]; \>> sort (A) \>> sort (A,2 ,'descend' ) \>> [X,I]=sort (A)
5.2 多项式计算
多项式的表示:n+1长度的行向量
多项式的四则预算:加减(相应向量加减)
P=Conv(p1,p2);相乘
[Q,r]=Deconv(p1,p2);除法;p1 = conv(Q,p2)+r
多项式的求导:polyder;
多项式的求值:polyval(p,x);ployvalm(p,x)
矩阵运算、点运算
多项式的求根:roots(p);给出全部根;poly函数由全部根建立多项式
例1 设f和g为两个多项式(见课件),求f(x)+g(x),f(x)-g(x),f(x)×g(x),f(x)/g(x)。
1 2 3 4 5 6 7 8 9 10 \>> f=[3 ,-5 ,0 ,-7 ,5 ,6 ]; g=[3 ,5 ,-3 ]; g1=[0 ,0 ,0 ,g]; \>> f+g1 \>> f-g1 \>> conv(f,g) \>> [Q,r]=deconv(f,g) \>> conv(g,Q)+r
例2 已知两个多项式a和b(见课件),计算两个多项式的乘积的导函数、商的导函数。
1 2 3 4 5 6 7 \>> a=[3 1 0 -6 ]; \>> b=[1 2 ]; \>> polyder(a) \>> c=polyder(a,b) \>> [p,q]=polyder(a,b)
>> [p,q]=polyder(a,b)
例3 以一个多项式为例,取一个2×2矩阵为自变量,分别用polyval和polyvalm计算该多项式的值。
1 2 3 4 5 \>> a=[1 ,8 ,0 ,0 ,-10 ]; \>> x=[-1 ,1.2 ;2 ,-1.8 ]; \>> y1=polyval(a,x) \>> y2=polyvalm(a,x)
例4 请推算空气流量在[0, 2]范围内什么水平时,加热效率为最高。
1 2 3 4 5 6 7 8 9 \>> p=[-38.89 ,126.11 ,-3.42 ]; \>> q=polyder(p) \>> roots(q) \>> polyval(p,1.6214 ) \>> x=0 :0.1 :2 ; \>> plot (x,polyval(p,x),1.6214 ,98.8154 ,'rp' );
5.3 数据插值(函数逼近的方法)
数值计算方法
Interp1(x,y,x1,method)
Method:linear(线性插值);nearest(最近点插值);pchip(分段3次埃尔米特插值);spline(3次样条插值);
Interp2(x,y,z,x1,y1,method)
引例 零件加工问题
1 2 3 4 5 \>> x= [0 ,3 ,5 ,7 ,9 ,11 ,12 ,13 ,14 ,15 ]; \>> y=[0 ,1.2 ,1.7 ,2.0 ,2.1 ,2.0 ,1.8 ,1.2 ,1.0 ,1.6 ]; \>> x1=0 :0.1 :15 ; \>> y1=interp1(x,y,x1,'spline' ); \>> plot (x1,y1)
粮储仓的通风控制问题
1 2 3 4 5 6 7 8 >> x=20 :10 :90 ; \>> y=(0 :5 :20 )'; \>> z=[8.9 ,10.32 ,11.3 ,12.5 ,13.9 ,15.3 ,17.8 ,21.3 ;8.7 ,10.8 ,11 ,12.1 ,13.2 ,14.8 ,16.55 ,20.8 ;8.3 ,9.65 ,10.88 ,12 ,13.2 ,14.6 ,16.4 ,20.5 ;8.1 ,9.4 ,10.7 ,11.9 ,13.1 ,14.5 ,16.2 ,20.3 ;8.1 ,9.2 ,10.8 ,12 ,13.2 ,14.8 ,16.9 ,20.9 ]; \>> xi=20 :90 ; \>> yi=(0 :20 )'; \>> zi=interp2(x,y,z,xi,yi,'spline' ); \>> surf(xi,yi,zi)
5.4 数据插值应用举例
机动车刹车距离问题
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 v=20 :10 :150 ; vs=v.*(1000 /3600 ); d1=10. *vs; d2=[3.15 ,7.08 ,12.59 ,19.68 ,28.34 ,38.57 ,50.4 ,63.75 ,78.71 ,95.22 ,113.29 ,132.93 ,154.12 ,176.87 ]; d3=10 ; d=d1+d2+d3; vi=20 :1 :150 ; di=interp1(v,d,vi,'spline' ); x=abs (di-120 ); [y,i ]=sort (x); vi(i (1 )) plot (vi,di,vi(i (1 )),di(i (1 )),'rp' ) \>> j =find (vi==125 ); \>> di(j ) \>> plot (vi,di,125 ,480.14 ,'rp' )
沙盘制作问题
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 x=0 :200 :1800 ; y=x'; z=[2000 ,2000 ,2001 ,1992 ,1954 ,1938 ,1972 ,1995 ,1999 ,1999 ;2000 ,2002 ,2006 ,1908 ,1533 ,1381 ,1728 ,1959 ,1998 ,2000 ; 2000 ,2005 ,2043 ,1921 ,977 ,897 ,1310 ,1930 ,2003 ,2000 ;1997 ,1978 ,2009 ,2463 ,2374 ,1445 ,1931 ,2209 ,2050 ,2003 ; 1992 ,1892 ,1566 ,1971 ,2768 ,2111 ,2653 ,2610 ,2121 ,2007 ;1991 ,1875 ,1511 ,1556 ,2221 ,1986 ,2660 ,2601 ,2119 ,2007 ; 1996 ,1950 ,1797 ,2057 ,2849 ,2798 ,2608 ,2303 ,2052 ,2003 ;1999 ,1999 ,2079 ,2685 ,3390 ,3384 ,2781 ,2165 ,2016 ,2000 ; 2000 ,2002 ,2043 ,2271 ,2668 ,2668 ,2277 ,2049 ,2003 ,2000 ;2000 ,2000 ,2004 ,2027 ,2067 ,2067 ,2027 ,2004 ,2000 ,2000 ]; surf(x,y,z); x1=0 :100 :1800 ; y1=x1'; z1=interp2(x,y,z,x1,y1,'spline' ); surf(x1,y1,z1); x2=0 :50 :1800 ; y2=x2'; z2=interp2(x1,y1,z1,x2,y2,'spline' ); surf(x2,y2,z2); contour(x2,y2,z2,12 )
5.5 曲线拟合(一种函数逼近的方法)
Polyfit:求得最小二乘拟合多项式系数
原理:函数逼近;多项式函数;误差最小(最小二乘法/最小平方法);
曲线拟合的实现方法:问题分析;
(若前面系数为0,则表明拟合次数太高,需要纠正)
引例—人口增长预测问题
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 x=1790 :10 :2010 ; y=[3.9 ,5.3 ,7.2 ,9.6 ,12.9 ,17.1 ,23.2 ,31.4 ,38.6 , 50.2 ,63.0 ,76.0 ,92.0 ,105.7 ,122.8 ,131.7 ,150.7 ,179.3 ,203.2 ,226.5 ,248.7 ,281.4 ,308.7 ]; plot (x,y,'*' ); p=polyfit(x,y,3 ); polyval(p,2020 ) plot (x,y,'*' ,x,polyval(p,x)); \>> polyval(p,2016 ) \>> x=1950 :10 :2010 ; \>> y=[150.7 ,179.3 ,203.2 ,226.5 ,248.7 ,281.4 ,308.7 ]; \>> p=polyfit(x,y,3 ) \>> p=polyfit(x,y,2 ); \>> plot (x,y,'*' ,x,polyval(p,x)) \>> polyval(p,2016 ) \>> polyval(p,2020 )
>> x=1950:10:2010; >> y=[150.7,179.3,203.2,226.5,248.7,281.4,308.7]; >> p=polyfit(x,y,3)
>> p=polyfit(x,y,2); >> plot(x,y,’*’,x,polyval(p,x)) >> polyval(p,2016)
>> polyval(p,2020)
家庭储蓄规律问题
1 2 3 4 5 \>> x=[0.6 ,1.0 ,1.4 ,1.8 ,2.2 ,2.6 ,3.0 ,3.4 ,3.8 ,4 ]; \>> y=[0.08 ,0.22 ,0.31 ,0.4 ,0.48 ,0.56 ,0.67 ,0.75 ,0.8 ,1.0 ]; \>> p=polyfit(x,y,1 ) \>> plot (x,y,'*' ,x,polyval(p,x))
5.6 曲线拟合应用举例 股票预测问题
1 2 3 4 5 6 7 8 9 10 \>> x=[2 ,3 ,4 ,5 ,8 ,9 ,10 ,11 ,12 ,15 ,16 ,17 ,18 ,19 ,22 ,23 ,24 ,25 ,26 ,29 ,30 ]; \>> y=[7.74 ,7.84 ,7.82 ,7.78 ,7.91 ,7.97 ,7.9 ,7.76 ,7.9 ,8.04 ,8.06 ,8.11 ,8.08 ,8.13 ,8.03 ,8.01 ,8.06 ,8.0 ,8.3 ,8.41 ,8.28 ]; \>> plot (x,y, '*' ); \>> p=polyfit(x,y,3 ); \>> plot (x,y,'*' ,x,polyval(p,x)); \>> x1=[31 ,32 ,33 ]; \>> xi=[x,x1]; \>> plot (x,y,'*' ,xi,polyval(p,xi)); \>> y1=[8.27 ,8.17 ,9.54 ]; \>> plot (x,y,'*' ,xi,polyval(p,xi),x1,y1,'rp' );
算法的参数优化问题
1 2 3 4 5 6 7 8 9 10 11 12 13 x=0.03 :0.03 :0.3 ; y1=[0.01 ,0.01 ,0.02 ,0.03 ,0.06 ,0.07 ,0.13 ,0.17 ,0.25 ,0.37 ]; y2=[0.85 ,0.76 ,0.68 ,0.62 ,0.56 ,0.52 ,0.49 ,0.46 ,0.43 ,0.39 ]; plot (x,y1,'*' ,x,y2,'o' );p1=polyfit(x,y1,2 ); p2=polyfit(x,y2,2 ); p=p1-p2; xi=roots(p); xj=0 :0.03 :0.36 ; yj1=polyval(p1,xj); yj2=polyval(p2,xj); yi=polyval(p1,xi(2 )) plot (x,y1,'*' ,x,y2,'o' ,xj,yj1,xj,yj2,xi(2 ),yi,'rp' );
专题五总结 6数值微积分与方程求解 6.1数值微分与数值积分
数值微分:
数值差分:向前差分;向后差分;中心差分;;diff
差商:向前差商;向后差商;中心差商;;
数值积分:牛顿-莱布尼兹公式;梯形法;辛普森法;高斯求积公式;
例1 设f(x)=sinx,在[0,2π]范围内随机采样,计算f’(x)的近似值,并与理论值f’(x)=cosx 进行比较。
1 2 3 4 5 6 \>> x=[0 ,sort (2 *pi *rand (1 ,5000 )),2 *pi ]; \>> y=sin (x); \>> f1=diff(y)./diff(x); \>> f2=cos (x(1 :end -1 )); \>> plot (x(1 :end -1 ),f1,x(1 :end -1 ),f2); \>> d=norm(f1-f2)
例2 分别用quad函数和quadl函数求定积分的近似值,并在相同的积分精度下,比较被积函数的调用次数。
1 2 3 4 5 6 7 8 9 \>> format long \>> f=@(x) 4. /(1 +x.^2 ); \>> [I,n]=quad(f,0 ,1 ,1e-8 ) \>> [I,n]=quadl(f,0 ,1 ,1e-8 ) \>> (atan (1 )-atan (0 ))*4 \>> format short
例3 求定积分。 被积函数文件fe.m:
1 2 3 4 5 6 function f =fe (x) f=1. /(x.*sqrt (1 -log (x).^2 )); \>> I=integral(@fe,1 ,exp (1 ))
例4 求定积分。 被积函数文件fsx.m:
1 2 3 4 5 6 function f =fsx (x) f=sin (1. /x)./x.^2 ; \>> I=quadgk(@fsx,2 /pi ,+Inf)
例5 设x=1:6,y=[6,8,11,7,5,2],用trapz函数计算定积分。
1 2 3 4 5 6 7 8 \>> x=1 :6 ; \>> y=[6 ,8 ,11 ,7 ,5 ,2 ]; \>> plot (x,y,'-ko' ); \>> grid on \>> axis([1 ,6 ,0 ,11 ]); \>> I1=trapz(x,y) \>> I2=sum(diff(x).*(y(1 :end -1 )+y(2 :end ))/2 )
例6 分别求二重积分和三重积分。
1 2 3 4 5 \>> f1=@(x,y) exp (-x.^2 /2 ).*sin (x.^2 +y); \>> I1=quad2d(f1,-2 ,2 ,-1 ,1 ) \>> f2=@(x,y,z) 4 *x.*z.*exp (-z.*z.*y-x.*x); \>> I2=integral3(f2,0 ,pi ,0 ,pi ,0 ,1 )
6.2线性方程组求解 直接法:高斯消去法;列主消元法(左除);矩阵的三角分解法;lu分解
例1 用左除运算符求解下列线性方程组。
1 2 3 \>> A=[2 ,1 ,-5 ,1 ;1 ,-5 ,0 ,7 ;0 ,2 ,1 ,-1 ;1 ,6 ,-1 ,-4 ]; \>> b=[13 ,-9 ,6 ,0 ]'; \>> x=A\b
例2 用LU分解求解例1中的线性方程组。
1 2 3 4 \>> A=[2 ,1 ,-5 ,1 ;1 ,-5 ,0 ,7 ;0 ,2 ,1 ,-1 ;1 ,6 ,-1 ,-4 ]; \>> b=[13 ,-9 ,6 ,0 ]'; \>> [L,U]=lu(A); \>> x=U\(L\b)
雅可比迭代法的函数文件jacobi.m:
1 2 3 4 5 6 7 8 9 10 11 12 13 function [y,n] =jacobi (A,b,x0,ep) D=diag (diag (A)); L=-tril (A,-1 ); U=-triu (A,1 ); B=D\(L+U); f=D\b; y=B*x0+f; n=1 ; while norm(y-x0)>=ep x0=y; y=B*x0+f; n=n+1 ; end
Gauss-Serdel迭代法的函数文件gauseidel.m:
1 2 3 4 5 6 7 8 9 10 11 12 13 function [y,n] =gauseidel (A,b,x0,ep) D=diag (diag (A)); L=-tril (A,-1 ); U=-triu (A,1 ); B=(D-L)\U; f=(D-L)\b; y=B*x0+f; n=1 ; while norm(y-x0)>=ep x0=y; y=B*x0+f; n=n+1 ; end
例3 分别用雅可比迭代法和高斯-赛德尔迭代法求解线性方程组。设迭代初值为0,迭代精度为10^(-6)。
1 2 3 4 5 \>> A=[4 ,-2 ,-1 ;-2 ,4 ,3 ;-1 ,-3 ,3 ]; \>> b=[1 ,5 ,0 ]'; \>> [x,n]=jacobi(A,b,[0 ,0 ,0 ]',1.0e-6 ) \>> [x,n]=gauseidel(A,b,[0 ,0 ,0 ]',1.0e-6 )
例4 分别用雅可比迭代法和高斯-赛德尔迭代法求解下列线性方程组,看是否收敛。
1 2 3 4 5 \>> A=[1 ,2 ,-2 ;1 ,1 ,1 ;2 ,2 ,1 ]; \>> b=[9 ;7 ;6 ]; \>> [x,n]=jacobi(A,b,[0 ;0 ;0 ],1.0e-6 ) \>> [x,n]=gauseidel(A,b,[0 ;0 ;0 ],1.0e-6 )
6.3 线性方程组应用举例 1.平面桁架结构受力分析问题
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 alpha=sqrt (2 )/2 ; A=[0 ,1 ,0 ,0 ,0 ,-1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ; 0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ; alpha,0 ,0 ,-1 ,-alpha,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ; alpha,0 ,1 ,0 ,alpha,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ; 0 ,0 ,0 ,1 ,0 ,0 ,0 ,-1 ,0 ,0 ,0 ,0 ,0 ; 0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ; 0 ,0 ,0 ,0 ,alpha,1 ,0 ,0 ,-alpha,-1 ,0 ,0 ,0 ; 0 ,0 ,0 ,0 ,alpha,0 ,1 ,0 ,alpha,0 ,0 ,0 ,0 ; 0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ,-1 ; 0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,0 ,0 ; 0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,alpha,0 ,0 ,-alpha,0 ; 0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,alpha,0 ,1 ,alpha,0 ; 0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,alpha,1 ]; b=[0 ;10 ;0 ;0 ;0 ;0 ;0 ;15 ;0 ;20 ;0 ;0 ;0 ]; f=A\b; disp (f')
2.小行星运行轨道计算问题
1 2 3 4 5 6 7 8 xi=[1.02 ,0.87 ,0.67 ,0.44 ,0.16 ]; yi=[0.39 ,0.27 ,0.18 ,0.13 ,0.13 ]; A=zeros (length (xi)); for i =1 :length (xi) A(i ,:)=[xi(i )*xi(i ),2 *xi(i )*yi(i ),yi(i )*yi(i ),2 *xi(i ),2 *yi(i )]; end b=-ones (length (xi),1 ); ai=A\b
或者:
1 2 3 4 5 6 7 xi=[1.02 ,0.87 ,0.67 ,0.44 ,0.16 ]'; yi=[0.39 ,0.27 ,0.18 ,0.13 ,0.13 ]'; A=[xi.*xi,2 *xi.*yi,yi.*yi,2 *xi,2 *yi]; b=-ones (length (xi),1 ); ai=A\b \>> f=@(x,y) 2.4645 *x.^2 -0.8846 *x.*y+6.4917 *y.^2 -1.3638 *x-7.2016 *y+1 ; \>> h=ezplot(f,[-0.5 ,1.2 ,0 ,1.2 ]);
6.4 非线性方程求解与函数极值计算
代法;二分法;割线法
Fzero
Fsolve
求极小值
Fminbnd
Fminsearch
Fminunc
fmincon
例1 求f(x)=0在x0=-5和x0=1作为迭代初值时的根。
1 2 3 4 5 6 7 8 9 10 11 12 13 \>> f=@(x) x-1. /x+5 ; \>> x1=fzero(f,-5 ) \>> x2=fzero(f,1 ) \>> x3=fzero(f,0.1 ) \>> f=@(x) x-1. /x+5 ; \>> x1=fsolve(f,-5 ,optimset('Display' ,'off' )) \>> x2=fsolve(f,1 ,optimset('Display' ,'off' )) \>> x3=fsolve(f,0.1 ,optimset('Display' ,'off' ))
f(x)=x^2-1=0的根
1 2 3 4 5 6 7 8 9 10 f=@(x) x.^2 -1 ; x=[]; x0=-0.25 :0.001 :0.25 ; for x00=x0 x=[x,fzero(f,x00)]; end plot (x0,x,'-o' )xlabel('初值' ); ylabel('方程的根' ); axis([-0.25 ,0.25 ,-1 ,1 ])
例2 求下列方程组在(1,1,1)附近的解并对结果进行验证。
1 2 3 4 5 6 \>> f=@(x) [sin (x(1 ))+x(2 )+x(3 )^2 *exp (x(1 )),x(1 )+x(2 )+x(3 ),x(1 )*x(2 )*x(3 )]; \>> f([1 ,1 ,1 ]) \>> x=fsolve(f,[1 ,1 ,1 ],optimset('Display' ,'off' )) \>> f(x)
例3 求函数f(x)在区间(-10,-1)和(1,10)上的最小值点。
1 2 3 4 \>> f=@(x) x-1. /x+5 ; \>> [xmin,fmin]=fminbnd(f,-10 ,-1 ) \>> [xmin,fmin]=fminbnd(f,1 ,10 )
例4 求解有约束最优化问题。
1 2 3 4 5 6 7 f=@(x) 0.4 *x(2 )+x(1 )^2 +x(2 )^2 -x(1 )*x(2 )+1 /30 *x(1 )^3 ; x0=[0.5 ;0.5 ]; A=[-1 ,-0.5 ;-0.5 ,-1 ]; b=[-0.4 ;-0.5 ]; lb=[0 ;0 ]; option=optimset('Display' ,'off' ); [xmin,fmin]=fmincon(f,x0,A,b,[],[],lb,[],[],option)
例5 要使每周送货车的里程最小,仓库应建在xy平面的什么位置?
1 2 3 4 5 a=[10 ,30 ,16.667 ,0.555 ,22.2221 ]; b=[10 ,50 ,29 ,29.888 ,49.988 ]; c=[10 ,18 ,20 ,14 ,25 ]; f=@(x) sum(c.*sqrt ((x(1 )-a).^2 +(x(2 )-b).^2 )); [xmin,fmin]=fminsearch(f,[15 ,30 ])
例6 在例5中,如果由于地域的限制,仓库必须建在曲线y=x^2上,则它应该建在何处?
1 2 3 4 5 6 7 8 9 10 11 12 非线性约束的函数文件funny.m: function [c,h] =funny (x) c=[]; h=x(2 )-x(1 )^2 ; \>> a=[10 ,30 ,16.667 ,0.555 ,22.2221 ]; \>> b=[10 ,50 ,29 ,29.888 ,49.988 ]; \>> c=[10 ,18 ,20 ,14 ,25 ]; \>> f=@(x) sum(c.*sqrt ((x(1 )-a).^2 +(x(2 )-b).^2 )); \>> [xmin,fmin]=fmincon(f,[15 ,30 ],[],[],[],[],[],[],'funny' )
6.5 常微分方程数值求解 单步法;
多步法
Solver
例1 求解微分方程初值问题,并与精确解进行比较。
1 2 3 4 f=@(t,y) (y^2 -t-2 )/4 /(t+1 ); [t,y]=ode23(f,[0 ,10 ],2 ); y1=sqrt (t+1 )+1 ; plot (t,y,'b:' ,t,y1,'r' );
例2 已知一个二阶线性系统的微分方程,取a=2,绘制系统的时间响应曲线和相平面图。
1 2 3 4 f=@(t,x) [-2 ,0 ;0 ,1 ]*[x(2 );x(1 )]; [t,x]=ode45(f,[0 ,20 ],[1 ,0 ]); subplot(2 ,2 ,1 );plot (t,x(:,2 )); subplot(2 ,2 ,2 );plot (x(:,2 ),x(:,1 ));
例3 假如点燃一个火柴,已知火焰球简化模型,分析λ的大小对方程求解过程的影响。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 lamda=0.01 ; f=@(t,y) y^2 -y^3 ; tic;[t,y]=ode45(f,[0 ,2 /lamda],lamda);toc disp ([‘ode45计算的点数' num2str(length(t))]); lamda=1e-5; f=@(t,y) y^2-y^3; tic;[t,y]=ode45(f,[0,2/lamda],lamda);toc disp([' ode45计算的点数' num2str(length(t))]); lamda=1e-5; f=@(t,y) y^2-y^3; tic;[t,y]=ode15s(f,[0,2/lamda],lamda);toc disp([' ode15s计算的点数' num2str(length(t))]);
6.6 常微分方程应用举例 1.Lotka-Volterra模型 第①问:
1 2 3 4 5 6 7 8 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,30 ],[300 ,150 ]) subplot(1 ,2 ,1 );plot (t,x(:,1 ),'-' ,t,x(:,2 ),'-*' ); legend ('x1(t)' ,'x2(t)' );xlabel('时间' );ylabel('物种数量' ); grid on subplot(1 ,2 ,2 );plot (x(:,1 ),x(:,2 )) grid on
第④问: 取λ=0.01, 所以稳定平衡点(1/λ,2/λ)即是(100,200),以此点作为初值时,画出其图像。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,30 ],[100 ,200 ]); plot (t,x(:,1 ),'-o' ,t,x(:,2 ),'-*' );legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' );
当将初始值变为(98,195)时,即向下十分接近平衡点,画出其图像。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,30 ],[98 ,195 ]); plot (t,x(:,1 ),'-o' ,t,x(:,2 ),'-*' );legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' );
当将初始值变为(70,150)时(向下偏离平衡点比较远时),画出其图像。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,30 ],[70 ,150 ]); plot (t,x(:,1 ),'-o' ,t,x(:,2 ),'-*' );legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' );
当将初始值变为(900,1600)时(向上偏离平衡点十分远时),画出其图像。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,500 ],[900 ,1600 ]); plot (t,x(:,1 ),t,x(:,2 ));legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' );
Lotka-Volterra改进模型 第①问:在原模型下,绘制狐狸和兔子数量的时间函数曲线。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; plot (t,x(:,1 ),t,x(:,2 ));legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' ); title('原模型下,狐狸和兔子数量的函数曲线' );
第②问:在改进模型下,狐狸和兔子数量的时间函数曲线。
1 2 3 4 5 6 7 8 rabbitFox=@(t,x) [2 *x(1 )*(1 -x(1 )/400 -0.005 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,50 ],[300 ,150 ]); plot (t,x(:,1 ),t,x(:,2 ));legend ('x1(t)-兔子' ,'x2(t)-狐狸' );xlabel('时间' ); ylabel('物种数量' ); title('改进模型下,狐狸和兔子数量的函数曲线' );
第③问:在原模型下,绘制狐狸数量相对于兔子数量的关系曲线。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [x(1 )*(2 -0.01 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,50 ],[300 ,150 ]); plot (x(:,1 ),x(:,2 ));xlabel('兔子数量' ); ylabel('狐狸数量' ); title('原模型下,狐狸数量相对于兔子数量的关系曲线' );
第④问:在改进模型下,狐狸数量相对于兔子数量的关系曲线。
1 2 3 4 5 6 7 rabbitFox=@(t,x) [2 *x(1 )*(1 -x(1 )/400 -0.005 *x(2 ));... x(2 )*(-1 +0.01 *x(1 ))]; [t,x]=ode45(rabbitFox,[0 ,50 ],[300 ,150 ]); plot (x(:,1 ),x(:,2 ));xlabel('兔子数量' ); ylabel('狐狸数量' ); title('改进模型下,狐狸数量相对于兔子数量的关系曲线' );
专题六总结 7 matlab 符号计算 7.1 符号对象: symbol
符号计算的结果是一个精确的数学表达式:符号推演,精确解。
数值计算的结果是一个数值:有舍入误差,近似解。
eval(): 将符号表达式转化成数值结果。
syms命令,定义符号变量
assume(condition)
assume(expr,set)
符号对象的运算:
四则运算
关系运算
逻辑运算
因式分解和展开运算:factor(s);expand(s);collect(s);collect(s,v)
其他运算
符号矩阵
梅森素数的验证问题。请验证M19 、 M23 、 M29 、 M31是否为梅森素数。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 >> syms p; \>> m=2 ^p-1 ; \>> p=19 ; \>> m19=eval(m) \>> factor (m19) \>> p=23 ; \>> m23=eval(m) \>> factor (m23) \>> p=29 ; \>> m29=eval(m) \>> factor (m29) \>> p=31 ; \>> m31=eval(m) \>> factor (m31)
例1 求方程ax^2+bx+c=0的根。
1 2 3 4 5 >> syms a b c x; \>> f=a*x^2 +b*x+c \>> g=coeffs(f,x) \>> g=g(end :-1 :1 ) \>> roots(g)
例2 当λ取何值时,以下齐次线性方程组有非零解。
1 2 3 4 syms lamda; A=[1 -lamda,-2 ,4 ;2 ,3 -lamda,1 ;1 ,1 ,1 -lamda]; D=det(A); factor (D)
7.2 符号微积分
极限:limit(f,x,a,’right/left’)
倒数:diff(f,x,n)
积分:
不定积分:int(f,x)
定积分:int(f,x,a,b); 无限:inf
例1 求下列极限。
1 2 3 4 5 >> syms a m x n; \>> f=(x^(1 /m)-a^(1 /m))/(x-a); \>> limit(f,x,a) \>> g=(1 +1 /n)^n; \>> limit(g,n,inf )
例2 求下列函数的导数。
1 2 3 4 5 6 7 8 >> syms x y z; \>> f=sqrt (1 +exp (x)); \>> diff(f) \>> g=x*exp (y)/y^2 ; \>> diff(g,x) \>> diff(g,y)
例3 求下列不定积分。
1 2 3 4 5 >> syms x t; \>> f=(3 -x^2 )^3 ; \>> int(f) \>> g=5 *x*t/(1 +x^2 ); \>> int(g,t)
例4 求下列定积分。
1 2 3 4 5 6 >> syms x t; \>> int(abs (1 -x),1 ,2 ) \>> int(1 /(1 +x^2 ),-inf ,inf ) \>> int(4 *x/t,t,2 ,sin (x))
河道水流量的估算问题。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 xi=0 :50 :600 ; yi=[4.4 ,4.5 ,4.6 ,4.8 ,4.9 ,5.1 ,5.4 ,5.2 ,5.5 ,5.2 ,4.9 ,4.8 ,4.7 ]; p=polyfit(xi,yi,3 ); plot (xi,yi,'o' ,xi,polyval(p,xi));syms y x; y=poly2sym(p,x); s=int(y,x,0 ,600 ); v=s*0.6 ; eval(v) yn=-yi; p=polyfit(xi,yn,3 ); plot (xi,yn,'o' ,xi,polyval(p,xi));syms y x yii; y=poly2sym(p,x); yii=diff(y,x); x=50 :60 ; k=eval(yii); all(abs (k)<1 /1.5 )
7.3 级数
级数求和
泰勒级数:将任意函数表示为幂级数
复杂函数的计算方法问题
加法= 逻辑电路
减法 = 加法转换
乘除 = 加法 + 移位
复杂函数 = 转换成四则运算,例泰勒级数
例1 求下列级数之和。
1 2 3 4 5 6 >> syms n; \>> s1=symsum(n^2 ,1 ,100 ) \>> s2=symsum((-1 )^(n-1 )/n,1 ,inf ) \>> s3=symsum((-1 )^(n-1 )/(2 *n-1 ),n,1 ,inf ) \>> hypergeom([-1 /2 , 1 ], 1 /2 , -1 ) \>> eval(s3)*4
银行利率的计算问题。
1 2 3 4 5 6 7 8 9 >> syms k r; \>> p2=symsum(50000 *(1 +0.045 /k)^k,k,2 ,2 ); \>> eval(p2) \>> p4=symsum(50000 *(1 +0.045 /k)^k,k,4 ,4 ); \>> eval(p4) \>> p12=symsum(50000 *(1 +0.045 /k)^k,k,12 ,12 ); \>> eval(p12) \>> limit((1 +r/k)^k,k,inf ) \>> 50000 *exp (0.045 )
例2 求函数f(x)在x=1处的5阶泰勒级数展开式。
1 2 3 4 >> syms x; \>> f=(1 +x+x^2 )/(1 -x+x^2 ); \>> taylor(f,x,1 ,'Order' ,6 ) \>>expand(ans )
例3 利用泰勒展开式计算三角函数的值。
1 2 3 4 5 >> syms x; \>> f=taylor(cos (x),x,pi ) \>> x=3 ; \>> eval(f) \>> cos (3 )
7.4 符号方程求解
代数方程:
solve(s) ;不一定准确
solve(s,v)
solve(s1,sn,v1,vn)
常微分方程 equation;condition;variable
D表示倒数: Dy;D2y
dsolve(e,c,v)
dsolve(e1,e2,c1,cn,v)
例1 解方程ax^2+bx+c=0。
1 2 3 4 5 6 7 >> syms x y a b c; \>> solve(a*x^2 +b*x+c==0 ) \>> f=a*x^2 +b*x+c==0 ; \>> solve(f) \>> solve(a*x^2 +b*x+c) \>> f=a*x^2 +b*x+c \>> solve(f)
例2 求下列微分方程或方程组的解。
1 2 3 >> syms x y t; \>> y=dsolve('Dy-(x^2+y^2)/x^2/2' ,x) \>> [x,y]=dsolve('Dx=4*x-2*y' ,'Dy=2*x-y' ,t)
疾病传染问题。
1 2 3 >> syms a b c y t; \>> f=dsolve('Dy=a*y*(1-y)-c*y' , 'y(0)=b' ,t) \>> f=dsolve('Dy=a*y*(1-y)-a*y' , 'y(0)=b' ,t)
专题七知识点总结 8 matlab 图形用户界面设计 8.1 图形窗口与坐标轴
图形对象:曲线、曲面、坐标轴、图形窗口
图形窗口对象
坐标轴对象
例1 绘制多个图形,并保存图形句柄。
1 2 3 4 5 6 7 8 t=0 :pi /10 :2 *pi ; h1=plot3 (t+pi ,t-2 *pi ,sin (t),'r' ); hold on;[x,y]=meshgrid (t); z=sin (x); h2=mesh(t-2 *pi ,t+pi ,z); [x3,y3,z3]=cylinder(t); h3=surf(x3,y3,z3);
例2 分别在两个子图中绘制曲线和曲面。然后设置子图1的背景色为黄色,曲线线条颜色为红色,设置子图2的背景色为青色。
1 2 3 4 5 6 7 8 9 10 subplot(1 ,2 ,1 ) h1=fplot(@(t)t.*sin (t),@(t)t.*cos (t),[0 ,6 *pi ] ); axis equal subplot(1 ,2 ,2 ) [x,y,z]=peaks(20 ); h2=mesh(x,y,z); h10=h1.Parent; h10.Color='y' ; h1.Color='r' ; h2.Parent.Color='cyan' ;
例3 建立一个图形窗口。该图形窗口没有菜单条,标题名称为“图形窗口示例”,起始于屏幕左下角、宽度和高度分别为300像素点和150像素点,背景颜色为青色,且当用户从键盘按下任意一个键时,然后在窗口中单击鼠标左键,在鼠标指针所在位置将显示“Hello,World!” 。
1 2 3 4 5 6 7 hf=figure ; hf.Color=[0 ,1 ,1 ]; hf.Position=[1 ,1 ,300 ,150 ]; hf.Name='图形窗口示例' ; hf.NumberTitle='off' ; hf.MenuBar='none' ; hf.ButtonDownFcn='gtext(''Hello,World!'')' ;
例4 利用坐标轴对象实现图形窗口的分割。
1 2 3 4 5 6 7 8 ha1=axes('Position' ,[0.1 ,0.1 ,0.7 ,0.7 ]); contour(peaks(20 )) ha1.Title=title('等高线' ); ha1.YLabel=ylabel('南北向' ); ha1.XLabel=xlabel('东西向' ); ha2=axes('Position' ,[0.65 ,0.7 ,0.28 ,0.28 ]); surf(peaks(20 )) ha2.View=[-30 ,45 ];
例5 定义包含4种颜色的ColorOrder属性,绘制6条曲线。
1 2 3 4 5 6 7 x=[0 ,0 ];y=[0 ,1 ]; ha=axes; ha.ColorOrder=[0 ,0 ,0 ; 1 ,0 ,0 ; 0 ,1 ,0 ; 0 ,0 ,1 ]; hold onplot (x,y, x+0.5 ,y, x+1 ,y, x+1.5 ,y, x+2 ,y, x+2.5 ,y);ha.XLim=[-0.2 ,3 ]; ha.YLim=[-0.2 ,1.2 ];
8.2 曲线与曲面对象
曲线对象
line(x,y,z)
plot 函数每调用一次,就会刷新坐标轴,清空原有图形再绘制新的曲线
曲面对象
surface(x,y,z,c)
surf函数每调用一次,就会刷新坐标轴,清空原有图形再绘制新的曲面
关照处理
图形对象的反射特性
例1 利用曲线对象绘制五环图案。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 t=-0.1 : 0.1 : 2 *pi ; x=cos (t); y=sin (t); line(x,y,'Color' ,'b' ) line(x+1.2 ,y-1 ,'Color' ,'y' ) line(x+2.4 ,y,'Color' ,'k' ) line(x+3.6 ,y-1 ,'Color' ,'g' ) line(x+4.8 ,y,'Color' ,'r' ) ha=gca; for n=1 :size (ha.Children) ha.Children(n).LineWidth=5 ; end ha.XLim=[-2 ,7 ]; ha.YLim=[-3 ,2 ]; axis equal
例2 利用曲面对象绘制立体圆环
1 2 3 4 5 6 7 8 9 10 11 r=linspace (0 ,2 *pi ,60 ); [u,v]=meshgrid (r); x=(8 +3 *cos (v)).*cos (u); y=(8 +3 *cos (v)).*sin (u); z=3 *sin (v); axes('view' ,[-37.5 ,30 ]) hs=surface(x,y,z); axis equal; hs.EdgeColor='none' ; hs.FaceColor='interp' ;
例3 绘制光照处理后的圆环面并观察不同光照模式下的效果。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 r=linspace (0 ,2 *pi ,60 ); [u,v]=meshgrid (r); x=(8 +3 *cos (v)).*cos (u); y=(8 +3 *cos (v)).*sin (u); z=3 *sin (v); axes('Position' ,[0.05 ,0.675 ,1.0 ,0.3 ], 'View' ,[-37.5 ,30 ]); hs1=surface(x,y,z); axis equal; hs1.EdgeColor='none' ; hs1.FaceColor='interp' ; axes('Position' ,[0.05 ,0.35 ,1.0 ,0.3 ], 'View' ,[-37.5 ,30 ]); hs2=surface(x,y,z); axis equal; hs2.EdgeColor='none' ; hs2.FaceColor='interp' ; light('Position' ,[0 ,0 ,8 ]) lighting flat axes('Position' ,[0.05 ,0.025 ,1.0 ,0.3 ], 'View' ,[-37.5 ,30 ]); hs3=surface(x,y,z); axis equal; hs3.EdgeColor='none' ; hs3.FaceColor='interp' ; light('Position' ,[0 ,0 ,8 ]) lighting phong
例4 绘制具有不同反射特性的圆环面并观察反射特性对图形效果的影响。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 r=linspace (0 ,2 *pi ,60 ); [u,v]=meshgrid (r); x=(8 +3 *cos (v)).*cos (u); y=(8 +3 *cos (v)).*sin (u); z=3 *sin (v); axes('Position' ,[0.05 ,0.675 ,1.0 ,0.3 ],'View' ,[-37.5 ,30 ]); hs1=surface(x,y,z);axis equal; hs1.EdgeColor='none' ;hs1.FaceColor='interp' ; light('Position' ,[0 ,0 ,8 ]) ;lighting phong hs1.SpecularStrength=0.1 ; axes('Position' , [0.05 ,0.35 ,1.0 ,0.3 ],'View' ,[-37.5 ,30 ]); hs2=surface(x,y,z);axis equal; hs2.EdgeColor='none' ;hs2.FaceColor='interp' ; light('Position' ,[0 ,0 ,8 ]) ;lighting phong hs2.SpecularStrength=0.5 ; axes('Position' , [0.05 ,0.025 ,1.0 ,0.3 ],'View' ,[-37.5 ,30 ]); hs3=surface(x,y,z);axis equal; hs3.EdgeColor='none' ;hs3.FaceColor='interp' ; light('Position' ,[0 ,0 ,8 ]) ;lighting phong ; hs3.SpecularStrength=1.0 ;
8.3 图形用户界面设计方法
Graphical User Interface(GUI)
控件函数
GUIDE
例1 在图形窗口中建立三个按钮对象,当单击按钮时分别绘制正弦曲线、显示或隐藏坐标轴的网格、清除坐标轴的图形。
1 2 3 4 5 6 7 8 9 10 11 ha= axes('Units' ,'pixels' ,'Position' ,[40 ,40 ,360 ,360 ]); ptgrid=uicontrol('Style' ,'pushbutton' ,... 'String' ,'网格' , 'Position' , [450 ,120 ,48 ,20 ],... 'Callback' ,'grid' ); btncla = uicontrol('Style' , 'pushbutton' , ... 'String' , '清除' ,'Position' , [450 ,80 ,48 ,20 ],... 'Callback' ,'cla' ); btnplot = uicontrol('Style' , 'pushbutton' , ... 'String' ,'绘图' ,'Position' , [450 ,160 ,48 ,20 ]); btnplot.Callback=@plot_sin;
新建回调函数文件plot_sin.m
1 2 3 4 function plot_sin (source, callbackdata) t=-pi :pi /20 :pi ; plot (t,sin (t)); end
例2 在例1的界面中添加“图形选项”菜单项,其中包括一个二级菜单项“线型”,其下又有3个子菜单项,分别为“实线”、“虚线”、“双划线”。
1 2 3 4 5 6 7 8 hopt=uimenu(gcf,'Label' ,'图形选项' ,'Accelerator' ,'L' ); hLStyle=uimenu(hopt,'Label' ,'线型' ,'Tag' ,'LStyle' , 'Enable' ,'off' ); hL_Solid=uimenu(hLStyle,'Label' ,'实线' ,... 'Tag' ,'Solid' ,'Callback' , @MLine_Type); hL_Dotted=uimenu(hLStyle,'Label' ,'虚线' ,... 'Tag' ,'Dotted' ,'Callback' , @MLine_Type); hL_Dashed=uimenu(hLStyle,'Label' ,'双划线' ,... 'Tag' ,'Dashed' ,'Callback' , @MLine_Type);
新建回调函数文件MLine_Type.m
1 2 3 4 5 6 7 8 9 10 function MLine_Type (source,callbackdata) hline=findobj('Type' ,'line' ); if strcmp(source.Tag,'Solid' )==1 hline.LineStyle='-' ; elseif strcmp(source.Tag,'Dotted' )==1 hline.LineStyle=':' ; elseif strcmp(source.Tag,'Dashed' )==1 hline.LineStyle='--' ; end end
修改回调函数文件plot_sin.m
1 2 3 4 5 6 > function plot_sin (source, callbackdata) > t=-pi :pi /20 :pi ; > plot (t,sin (t)); > h1=findobj('Tag' ,'LStyle' ); > h1.Enable='On' ; > end
8.4 用户界面设计工具
回调函数Callback
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 function pushbutton1_Callback (hObject, eventdata, handles) A=eval(handles.editfz.String); f=eval(handles.editpl.String)/50 ; theta=eval(handles.editxj.String)/180 *pi ; x=linspace (0 ,2 *pi ,60 ); if handles.OpSin.Value==1 y=A*sin (f*x+theta); else y=A*cos (f*x+theta); end plot (x,y);handles.PStyle.Enable='On' ; function Solid_Callback (hObject, eventdata, handles) hline=findobj('Type' ,'line' ); hline.LineStyle='-' ; handles.Solid.Checked='On' ; handles.Dotted.Checked='Off' ; handles.Dashed.Checked='Off' ; function r_Callback (hObject, eventdata, handles) hline=findobj('Type' ,'line' ); hline.Color='r' ; handles.r.Checked='On' ; handles.g.Checked='Off' ; handles.b.Checked='Off' ;
8.5 APP 设计工具
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 function Draw_Plot (app) f=eval(app.DiscreteKnob.Value); theta=app.Knob.Value/180 *pi ; x=linspace (0 ,2 *pi ,60 ); if app.RadioButton.Value==1 y=sin (f*x+theta); else y=cos (f*x+theta); end plot (app.UIAxes,x,y,'LineWidth' ,0.2 );end function Clear_Axes (app) cla(app.UIAxes) end
8.6 图形用户界面应用举例
例1 利用GUIDE设计工具设计如图所示的用户界面。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 handles.peaks=peaks(34 ); handles.membrane=membrane; [x,y]=meshgrid (-8 :0.3 :8 ); r=sqrt (x.^2 +y.^2 ); sinc=sin (r)./(r+eps ); handles.sinc=sinc; handles.current_data=handles.sinc; colormap(spring); mesh(handles.current_data) surf(handles.current_data) contour3(handles.current_data) if hObject.Value==1 grid on hObject.String='隐藏网格' ; else grid off hObject.String='显示网格' ; end str=hObject.String; val=hObject.Value; switch strtrim(str{val})case 'Peaks' handles.current_data=handles.peaks; case 'Membrane' handles.current_data=handles.membrane; case 'Sinc' handles.current_data=handles.sinc; end guidata(hObject,handles) str=hObject.String; cm=hObject.Value; colormap( eval(str{cm}) ); el=eval(handles.edit_el.String); az=eval(handles.edit_az.String); view(az,el) switch eventdata.NewValue.Tag case 'rb_flat' shading flat; case 'rb_interp' shading interp; case 'rb_faceted' shading faceted; end
例2 生成一个用于观察周期信号波形叠加效果的程序模块。
1 2 3 4 5 6 7 8 9 10 amr=eval(app.ARadio.Value); theta=app.PhDiff.Value/180 *pi ; x=linspace (0 ,10 *pi ,300 ); y=sin (x)+sin (3 *x+theta)/amr; if strcmp(app.Noise.Value,'On' ) y=awgn(y,30 ); end app.WaveA.Value=max (y); plot (app.UIAxes,x,y,'LineWidth' ,0.1 );
专题八总结 Simulink 系统仿真 9.1 simulink 仿真基础
启动
“主页”选项卡:新建
“主页”选项卡:simulink
命令行窗口,simulink
模型
模块
模型存盘
模块参数
步骤
建立系统仿真模型
设置仿真参数
启动仿真并分析仿真结果
9.2 子系统的创建与封装 创建
新建
模块转换
封装
条件执行
9.3 s函数的设计与应用
S函数: 系统函数,system function
采用S函数实现y=kx+b。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 function [sys,x0,str,ts] =timekb (t,x,u,flag,k,b) switch flag case 0 [sys,x0,str,ts]=mdlInitializeSizes; case 3 sys=mdlOutputs(t,x,u,k,b); case {1 ,2 ,4 ,9 } sys=[]; otherwise error(num2str(flag)) end function [sys,x0,str,ts] =mdlInitializeSizes () sizes=simsizes; sizes.NumContStates=0 ; sizes.NumDiscStates=0 ; sizes.NumOutputs=1 ; sizes.NumInputs=1 ; sizes.DirFeedthrough=1 ; sizes.NumSampleTimes=1 ; sys=simsizes(sizes); x0=[]; str=[]; ts=[-1 ,0 ]; function sys =mdlOutputs (t,x,u,k,b) sys=k*u+b;
9.4 simulink仿真应用举例 仿真求解:
(1)求最大安全体重
1 2 3 4 5 6 7 8 9 10 11 for m=100 :-0.5 :20 [t,x,y_w]=sim('bengji' ,0 :0.01 :100 ); if min (y_w)>1.5 break ; end end disp (['最大安全体重是' ,num2str(m)])dis=min (y_w); disp (['最小的安全距离是' ,num2str(dis)])plot (t,y_w)grid on
(2)求最小弹性系数
1 2 3 4 5 6 7 8 9 10 11 12 m=65 ; for k=10 :0.1 :50 [t,x,y_w]=sim('bengji' ,0 :0.01 :100 ); if min (y_w)>1.5 break ; end end disp (['最小弹性系数k是' ,num2str(k)])dis=min (y_w); disp (['最小的安全距离是' ,num2str(dis)])plot (t,y_w)grid on
专题九总结 10 外部程序接口 10.1 在Excel中使用matlab Spreadsheet Link
10.2 matlab 文件操作
fid = fopen(filename,permission)
fid, 文件识别号,不成功时为-1
permission,’r’,只读;’w’,写;’a’,追加;‘r+’,读写。
status=fclose(fid)
fid,文件标识号,为all时,关闭已所有打开的文件
status,0,关闭成功;-1,关闭不成功
读写文本文件
[A,count]=fscanf(fid,fmt,size)
count=fprintf(fid,fmt,A)
读写二进制文件
[A,count]=fread(fid,size,precision,skip)
count=fwrite(fid,A,precious)
数据文件定位:文件位置指针
fseek(fid,offset,origin)
ftell(fid)
status=feof(fid)
例1 文件“观测记录.txt”的内容如图所示,读取文件前10行的数据。
1 2 3 4 5 6 7 8 9 10 11 12 fid=fopen('观测记录.txt' ,'r' ); title=fscanf(fid,'%s' ,6 ); qxsj=[]; for i =1 :10 qxsj{i ,1 }=fscanf(fid,'%s' ,1 ); qxsj{i ,2 }=fscanf(fid,'%s' ,1 ); qxsj{i ,3 }=fscanf(fid,'%f' ,1 ); qxsj{i ,4 }=fscanf(fid,'%f' ,1 ); qxsj{i ,5 }=fscanf(fid,'%f' ,1 ); qxsj{i ,6 }=fscanf(fid,'%s' ,1 ); end fclose(fid);
例2 计算y=exp(x)sinx,其中x∈[0,2π]。将x、y写入二进制文件“模拟数据.dat”。
1 2 3 4 5 fid=fopen('模拟数据.dat' ,'w' ); x=linspace (0 ,2 *pi ,100 ); y=exp (x).*sin (x); count=fwrite(fid, [x; y], 'double' ); fclose(fid);
例3 读取例2生成的文件中的后40组数据,并绘制图形。
1 2 3 4 5 6 7 8 9 fid=fopen('模拟数据.dat' ,'r' ); status=fseek(fid, -40 *2 *8 , 'eof' ); x=[]; y=[]; while ~feof(fid) x=[x; fread(fid,1 ,'double' )]; y=[y; fread(fid,1 ,'double' )]; end plot (x,y);fclose(fid);
10.3 在其他语言程序中读写matlab的数据文件
接口函数
mat文件 :matlab存储数据的标准格式;变量的名称、类型和值均保存。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 #include "mat.h" \#include <iostream> using namespace std;int main () { MATFile *pmat; mxArray *pa1, *pa2, *pa3; double data[9 ] = {1.1 , 2.2 , 3.3 , 4.4 , 5.5 , 6.6 , 7.7 , 8.8 , 9.9 }; const char *file = "mattest.mat" ; int status; cout << "生成文件 : " << file << endl; pmat = matOpen (file, "w" ); if (pmat == NULL ) { cout << "不能创建文件 : " << file << endl; cout << "(请确认是否有权限访问指定文件夹?)\n" ; return (0 ); } pa1 = mxCreateDoubleScalar (1.234 ); if (pa1 == NULL ) { cout << "不能创建变量。\n" ; return (0 ); } pa2 = mxCreateDoubleMatrix (3 , 3 , mxREAL); if (pa2 == NULL ) { cout << "不能创建矩阵。\n" ; return (0 ); } memcpy ((void *)(mxGetPr (pa2)), (void *)data, sizeof (data)); pa3 = mxCreateString ("MAT文件示例" ); if (pa3 == NULL ) { cout << "不能创建字符串。\n" ; return (0 ); } status = matPutVariable (pmat, "LocalDouble" , pa1); if (status != 0 ) { cout << "写入局部变量时发生错误。\n " ; return (0 ); } status = matPutVariableAsGlobal (pmat, "GlobalDouble" , pa2); if (status != 0 ) { cout << "写入全局变量时发生错误。\n" ; return (0 ); } status = matPutVariable (pmat, "LocalString" , pa3); if (status != 0 ) { cout << ": 写入String类型数据时发生错误。\n " ; return (0 ); } mxDestroyArray (pa1); mxDestroyArray (pa2); mxDestroyArray (pa3); if (matClose (pmat) != 0 ) { cout << "关闭文件时发生错误 " << endl; return (0 ); } cout << "文件创建成功!\n" ; return (1 ); }
10.4 在matlab中调用其他语言编写的程序 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 #include "mex.h" bool isCoprime (double *inx,double *outy) { int x=*inx,y=*outy,tmp; if (x<y) {tmp=x; x=y; y=tmp;} tmp=x%y; while (tmp) { x=y; y=tmp; tmp=x%y; } if (y==1 ) return true ; else return false ; } void mexFunction ( int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[] ) { double *result; int m,n,i; if (nrhs!=2 ) { mexErrMsgTxt ("输入参数应有两个!" ); return ; } if (nlhs!=1 ) { mexErrMsgTxt ("输出参数应只有1个!" ); return ; } for (i=0 ;i<2 ;i++){ m = int (mxGetM (prhs[i])); n = int (mxGetN (prhs[i])); if ( !mxIsDouble (prhs[i]) || m!=1 || n!=1 ) { mexErrMsgTxt ("输入参数必须是单个的数!" ); return ; } } plhs[0 ] = mxCreateDoubleMatrix (m,n, mxREAL); result = mxGetPr (plhs[0 ]); if (isCoprime (mxGetPr (prhs[0 ]),mxGetPr (prhs[1 ]))) *result=1 ; else *result=0 ; }
10.5 在其他语言程序中调用matlab函数 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 #include <engine.h> \#include <iostream> using namespace std;int main () { Engine *ep; mxArray *T = NULL , *R = NULL ; double t[180 ],r[180 ]; double a,b; a=2 ;b=3 ; for (int i=0 ;i<180 ;i++) { t[i]=i*0.1 ; r[i]=a+b*t[i]; } if (!(ep = engOpen ("" ))) { cout<< "\n不能启动MATLAB引擎\n" ; return 0 ; } T = mxCreateDoubleMatrix (1 , 180 , mxREAL); memcpy ((void *)mxGetPr (T), (void *)t, sizeof (t)); engPutVariable (ep, "T" , T); R = mxCreateDoubleMatrix (1 , 180 , mxREAL); memcpy ((void *)mxGetPr (R), (void *)r, sizeof (r)); engPutVariable (ep, "R" , R); engEvalString (ep, "polar(T,R);" ); engEvalString (ep, "title('阿基米德螺线');" ); system ("pause" ); mxDestroyArray (T); mxDestroyArray (R); engClose (ep); return 1 ; }
专题十总结 其他学习资源 https://matlabacademy.mathworks.com/cn
https://ww2.mathworks.cn/videos/getting-started-with-matlab-68985.html