MATLAB

《科学计算与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
p=[1,-3,1];
x=roots(p)

绘图:

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)-2i;
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; %递归调用求(n-1)!
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 BETA
f=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=[-1*z(2,:);z(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('物种数量');
  1. 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 符号微积分

  1. 极限:limit(f,x,a,’right/left’)

  2. 倒数:diff(f,x,n)

  3. 积分:

    • 不定积分: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. 级数求和
    • sum()
    • symsum(s,v,n,m)
  2. 泰勒级数:将任意函数表示为幂级数

    • taylor(f,v,a,Name,Value)
  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 符号方程求解

  1. 代数方程:
    • solve(s) ;不一定准确
    • solve(s,v)
    • solve(s1,sn,v1,vn)
  2. 常微分方程 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. 图形对象:曲线、曲面、坐标轴、图形窗口

    • 句柄

    • 属性

      points 磅

  2. 图形窗口对象

  3. 坐标轴对象

例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 on
plot(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 曲线与曲面对象

  1. 曲线对象

    • line(x,y,z)

      plot 函数每调用一次,就会刷新坐标轴,清空原有图形再绘制新的曲线

  2. 曲面对象

    • surface(x,y,z,c)

      surf函数每调用一次,就会刷新坐标轴,清空原有图形再绘制新的曲面

  3. 关照处理

    • light()
  4. 图形对象的反射特性

例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)

  1. 控件函数
    • uicontrol()
  2. 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]);
%设置“绘图”按钮的Callback属性值是plot_sin函数句柄。
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 用户界面设计工具

  1. 回调函数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
%“绘图”按钮的Callback函数
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); %获取离散旋钮对象的Value属性值
theta=app.Knob.Value/180*pi; %获取旋钮对象的Value属性值
x=linspace(0,2*pi,60);
if app.RadioButton.Value==1
y=sin(f*x+theta);
else
y=cos(f*x+theta);
end
% plot函数的第1个参数必须是坐标轴对象句柄
%plot函数默认的线条宽度为1,App程序的线条宽度必须小于0.4
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);

%Mseh_Callback函数
mesh(handles.current_data)

%Surf_Callback函数
surf(handles.current_data)

%Contour3_Callback函数
contour3(handles.current_data)

%切换按钮的Callback函数
if hObject.Value==1
grid on
hObject.String='隐藏网格';
else
grid off
hObject.String='显示网格';
end

%列表框的Callback函数
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)

%色图弹出式菜单的Callback函数
str=hObject.String;
cm=hObject.Value;
colormap( eval(str{cm}) );

%视点设置按钮的Callback函数
el=eval(handles.edit_el.String);
az=eval(handles.edit_az.String);
view(az,el)

%着色方式按钮组ChooseShading的SelectionChanged函数
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 系统仿真

启动

  1. “主页”选项卡:新建
  2. “主页”选项卡:simulink
  3. 命令行窗口,simulink

模型

模块

模型存盘

模块参数

步骤

  1. 建立系统仿真模型
  2. 设置仿真参数
  3. 启动仿真并分析仿真结果

9.2 子系统的创建与封装

创建

  1. 新建
  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=[]; %将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 文件操作

  1. fid = fopen(filename,permission)
    • fid, 文件识别号,不成功时为-1
    • permission,’r’,只读;’w’,写;’a’,追加;‘r+’,读写。
  2. status=fclose(fid)
    • fid,文件标识号,为all时,关闭已所有打开的文件
    • status,0,关闭成功;-1,关闭不成功
  3. 读写文本文件
    • [A,count]=fscanf(fid,fmt,size)
    • count=fprintf(fid,fmt,A)
  4. 读写二进制文件
    • [A,count]=fread(fid,size,precision,skip)
    • count=fwrite(fid,A,precious)
  5. 数据文件定位:文件位置指针
    • 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存储数据的标准格式;变量的名称、类型和值均保存。

  • save
  • load
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; /* 定义MAT文件指针*/
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; //定义存储函数返回值的变量
/* 打开一个MAT文件,如果不存在则创建一个MAT文件,如果打开失败,则返回 */
cout << "生成文件 : " << file << endl;
pmat = matOpen(file, "w");
if (pmat == NULL) {
cout << "不能创建文件 : " << file << endl;
cout << "(请确认是否有权限访问指定文件夹?)\n";
return(0);
}
/* 创建3个mxArray对象,其中pa1存储一个实数,pa2为3×3的双精度实型矩阵,*/
/* pa3存储字符串,如果创建失败则返回 */
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);
}
/* 向MAT文件中写数据,失败则返回 */
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);
/* 关闭MAT文件 */
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) //最大公约数为1,所以互质
return true;
else //最大公约数大于1,所以不互质
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; }
}
//在MATLAB工作区建立一个矩阵,用于存放计算结果,
//矩阵的大小与输入实参的大小相同。
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,用于操作引擎对象。
Engine *ep;
//定义mxArray 类型的指针,用于指向所调用MATLAB函数的输入对象。
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]; }
//启动MATLAB计算引擎
//Windows系统中,engOpen函数的参数为空字符串,指定在本机启动计算引擎。
if (!(ep = engOpen(""))) {
cout<< "\n不能启动MATLAB引擎\n";
return 0;
}
//建立MATLAB变量,调用MATLAB函数
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