04 程序设计基础

Download Report

Transcript 04 程序设计基础 

程序设计基础
一、M脚本文件及编辑器
二、结构化程序设计
三、顺序结构：以线性拟合为例
四、分支(选择)结构：if, switch
五、循环结构：for, while
六、错误控制结构：try-catch
一、M脚本文件及编辑器
• 是一串按用户意图排列而成matlab指令集合
• 文件的扩展名为.m，文件为纯文本文件（可用记
事本等纯文本编辑器进行编辑）
• 文件按Maltab变量名的规则进行命名，不可包含
中文等字符，否则无法在Maltab中进行调用（虽
然可以正常存盘）。
• 脚本文件运行后，所产生的变量都驻留在matlab
的 基本工作空间（base workspace）中。
1、创建新的m脚本文件
也可在command windows中，输入：
>> edit newfilename.m
2、打开已存在的m脚本文件
也可在command windows中，输入：
>> edit filename.m
3、脚本文件的编辑器 Editor
行尾加分号，则这一行的执行结果不会在command
windows中显示，可以加快脚本的执行速度。(显示
大量的运算结果是很费时间的)
4、m脚本文件存盘
保存脚本文件，文件名 必须按Maltab变量名的规
则进行命名，否则无法在Maltab中进行调用
5、运行m脚本文件
方法1：在编辑器（editor）中打开m文件，点击运行图标，
一般在编写、调试m脚本文件的时候使用。
方法2：在command windows窗口中直接输入m脚本文件的文件
名（文件扩展名.m不能输） ，然后回车，即可运行脚本文件。
方法3：在current directory窗口中选择m脚本文件，然后点击鼠
标右键，再弹出的快捷菜单中，选择 run 。
M脚本文件的运行路径 ( path )
若M脚本文件不是保存在Matlab的当前工作目录下，则第一次
运行文件时会出现如下的对话框，此时选择”change Directory”或
“Add to Path”都可以。Path还可用“File”菜单下的“Set
Path…”命令进行设置。
6、使用编辑器的 Cell Mode
执行整个文件中的所有语句
执行光标所在cell内的所有语句，完成后光标跳至下一个cell
执行光标所在cell内的所有语句
若cell mode工具条上出现的按钮不全，可以右击工具
条，选择自定义功能进行自定义
恢复默认设置
可自定义在工
具条上显示哪
些按钮
二、结构化程序设计
• Matlab支持结构化程序设计：
(1)结构化程序：由三种基本结构反复嵌套构成的程序
(2)基本思想：任何程序都可以用三种基本结构表示
• 三种基本结构：
(1)顺序结构
(2)分支(选择)结构
(3)循环结构
• 特点：
(1)优点：结构清晰，易读，提高程序设计质量和效率
(2)缺点：用户需求难以准确定义，适应变化的能力弱，开发周
期长
三、顺序结构：以线性拟合为例
顺序结构：顺序执行程序的各条语句
输入：已知条件
处理
输出：结果
线性拟合(Linear Fit)模型
Y与X具有统计
关系而且是线性
建立
回归模型
Yi=β0+β1Xi+εi
(i=1,2,···,n)
其中：
(X i,Yj)表示(X,Y)的第i个观测值
β0 ,β1为参数，β0+β1Xi为反映统计关系直线的分量
εi为反映在统计关系直线周围散布的随机分量
εi～N (0,σ2)， εi 服从正态分布
Yi=β0+β1Xi+εi
β0和β1均未知
根据样本数据
对β0和β1
进行估计
β0和β1的估计
值为b0和b1
建立一元线性回归方程
Yˆ  b0  b1 X
一般而言，所求的b0和b1应能使每个样本观测点
(Xi,Yi)与回归直线之间的偏差尽可能小。
一元线性回归方程
最小
二乘法
Y与X之间
为线性关系
The least square method
参考资料：
最小二乘法（百度百科）
最小二乘法（维基百科）
谷神星、高斯、最小二乘法
选出一条最能反
映Y与X之间关系
规律的直线
n
令 Q   Yi   b0  b1 X i  
2
i 1
Q达到最小值
b0和b1称为最小二乘估计量
n
Q
 2 Yi   b0  b1 X i    0
b0
i 1
n
Q
 2 Yi   b0  b1 X i   X i  0
b1
i 1
n
b1 
 X
i 1
i
 X Yi  Y 
n
 X
i 1
i
X
2
n

 X
i 1
n
 X
i 1
 X  Yi
i
i
X
2
微积分极值
的必要条件
, b0  Y  b1 X
确定系数(Coefficient Of Determination)
残差 residuals ：ei  Yi  Yˆi  Yi  b0  b1 X i
残差越小，各观测值聚集在回归直线周围的紧密程度就越
大，说明直线与观测值的拟合越好，定义确定系数(COD)为：

n
R2 
i 1
n
Yˆi  Y

 Y  Y 
i 1
i
0  R 1
2
 1
2

n
2
i 1
n
Yi  Yˆi

 Y  Y 
i 1
i
n
2
e
i
2
 1
2
i 1
n
 Y  Y 
i 1
2
i
一般情况下，R2的值越大，拟合得越
好。
线性拟合的相关系数
（Correlation coefficient）
n
r R 
(X
i 1
2
i
 X )(Yi  Y )
n
2
(
X

X
)
 i
i 1
1  r  r
r 与斜率 b1 取相同的符号
r = 1：
r = -1：
r = 0：
完全正相关
完全负相关
无线性关系
n
2
(
Y

Y
)
 i
i 1
例：测得铜导线在温度Ti下的电阻为Ri，编写一个M脚本文件，对
数据进行线性拟合，求出斜率、截距和相关系数。
n
T(℃)
R(Ω)
1
19.1
76.30
2
25.0
77.80
3
30.1
79.25
4
36.0
80.80
5
40.0
82.35
6
45.1
83.90
7
50.0
85.10
R  b0  bT
1
1、如何输入已知条件？
温度 T、电阻R的测量数据可以用matlab的一维
数值数组存放，便于后继的处理。
清除当前命令窗口中显示的内容，以便于观
察程序的输出结果。
clc
清除内存中的所有变量，避免混淆。
clear all
x = [19.1,25.0,30.1,36.0,40.0,45.1,50.0];
y = [76.30,77.80,79.25,80.80,82.35,83.90,85.10];
2、如何处理？
把线性拟合的数学语言“翻译”成matlab语言
测量值：X 1 ,
Y1 ,
, Xn
, Yn
数学模型：Yˆ  b0  b1 X
n
最佳参数：b1 
 X
i 1
n
 X
i 1
 X  Yi
i
i
X
2
n
拟合结果评价：r 
(X
i 1
i
 X )(Yi  Y )
n
 (Xi  X )
i 1
b0  Y  b1 X
，
n
2
2
(
Y

Y
)
 i
i 1
数学语言翻译成Matlab语言
求和函数：sum()
n
b1 
 X
i 1
n
 X
i 1
 X  Yi
i
i
X
2
数组相乘 .*
n
， b0  Y  b1 X ， r 
数组求幂 .^
(X
i 1
n
i
 X )(Yi  Y )
2
(
X

X
)
 i
i 1
n
2
(
Y

Y
)
 i
i 1
xm = mean(x)
ym = mean(y)
b1 = sum((x-xm).*y)/sum((x-xm).^2)
b0 = ym - b1*xm
r = sum((x-xm).*(y-ym))/sqrt(sum((x-xm).^2)*sum((y-ym).^2))
3、如何输出结果？
用字符串以较为人性化的方式输出。
disp(x)
num2str(x)
[s1,s2]
在命令窗口显示x的值，x可以是数值、
字符串等
将数值变量x转换为字符串
把两个字符串s1，s2合并为一个长字符
串
disp('本次线性拟合的结果为： ')
disp(['截距b0 = ',num2str(b0)])
disp(['斜率b1 = ',num2str(b1)])
disp(['相关系数r = ',num2str(r)])
如何更人性化的输出？
用plot()函数绘图：实验散点图、
拟合直线图、图形注释等。
实现线性拟合的matlab完整程序
clc
clear all
x = [19.1,25.0,30.1,36.0,40.0,45.1,50.0];
y = [76.30,77.80,79.25,80.80,82.35,83.90,85.10];
xm = mean(x);
ym = mean(y);
b1 = sum((x-xm).*y)/sum((x-xm).^2);
b0 = ym - b1*xm;
r = sum((x-xm).*(y-ym))/sqrt(sum((x-xm).^2)*sum((y-ym).^2));
disp('本次线性拟合的结果为： ')
disp(['截距b0 = ',num2str(b0)])
disp(['斜率b1 = ',num2str(b1)])
disp(['相关系数r = ',num2str(r)])
4、使用输入函数input()改进程序，提高实用性
v = input(‘str’)
在屏幕上显示字符串str定义的提示信息 ，并等待用户的键
盘输入。可以从键盘输入任意的matlab表达式，按回车键确
认。
(1)如直接回车，将空数组[ ]赋值给变量v。
(2)若输入10,将数值10赋值给变量v
(3)若输入[1,2,3,4]，将一维数值数组[1,2,3,4]赋值给变量v
(4)若输入’10’,将字符串’10’赋值给变量v 。
(5)若输入不合法，matlab会显示错误信息，并继续在屏幕上显
改进的线性拟合程序
clc
clear all
x = input('请输入进行线性拟合的x数组：');
y = input('请输入进行线性拟合的y数组： ');
xm = mean(x);
ym = mean(y);
b1 = sum((x-xm).*y)/sum((x-xm).^2);
b0 = ym - b1*xm;
r = sum((x-xm).*(y-ym))/sqrt(sum((x-xm).^2)*sum((y-ym).^2));
disp('本次线性拟合的结果为： ')
disp(['截距b0 = ',num2str(b0)])
disp(['斜率b1 = ',num2str(b1)])
disp(['相关系数r = ',num2str(r)])
5、用Matlab提供的拟合函数进行线性拟合
clc
clear all
x = [19.1,25.0,30.1,36.0,40.0,45.1,50.0];
y = [76.30,77.80,79.25,80.80,82.35,83.90,85.10];
p = polyfit(x,y,1)
c = corrcoef(x,y)
b1 = p(1)
b0 = p(2)
r = c(1,2)
p = polyfit(x,y,n)：多项式拟合函数，n=1
时，就是线性拟合，返回值p是行数组。
p = [b1,b0]
c = corrcoef(x,y)：计算向量x，y的相关
度，返回值是2×2矩阵。
相关系数 r = c(1,2) = c(2,1)
四、分支(选择)结构
选择结构也称为决策结构、分支结构或判断结构
1、 if - end 结构
2、 if - else - end 结构
3、 if - elseif - … - else - end 结构
4、 switch 构
1、 if - end 结构
if expression
statements
end
这是最简单、最常
用的选择结构。
expression 为条件表达式 ， statements 为要执行
的语句。 expression的值为逻辑真（非0）时，statements
才会执行。
clc
n = input('请输入一个整数：');
str = '你输入的不是正偶数。';
if rem(n,2) == 0
str = '你输入的是正偶数。';
end
disp(str)
分别输入4 和 -4
看看输出结果对不对。
if expression
statements
end
一般情况下，表达式 expression 都是标量，但也允许
为数组，此时只有数组元素都为逻辑真时，statements才被
执行。如果表达式为空数组，被认为是假。
clc
n = input('请输入一个整数数组：');
str = '你输入的数不全是正偶数。';
分别输入[4, -2]和[4,2]
if rem(n,2)==0&n>0
看看输出结果对不对。
str = '你输入的数全部都是正偶数。';
end
disp(str)
2、 if - else – end 结构
若要在 expression 为 True和 False 两种条件下执
行不同的操作，可以使用如下格式的 if - else - end 结构：
if expression
statements 1
else
statements 2
end
3、 if - elseif - … - end 结构
当需要根据多个条件执行多个不同的操作时，可以采
用下面的选择结构，matlab 将从上到下检测各个表达式，执
行与所遇到的第一个为 True 的表达式相对应的命令集，
然后退出选择结构。
if expression1
statements1
elseif expression2
statements2
elseif expression3
statements3
……
else
statements
end
注意：
elseif 中间没有空格！
4、 switch 结构
switch expression
case expr1
statement1
case expr2
statement2
……
otherwise
statement
end
当遇到 switch 结构时，
MATLAB 将 expression
的值依次与各个 case 指令后
面的检测值进行比较：
若比较结果为假，则取下
一个检测值再比较；
若比较结果为真，则执行
相应的一组指令，然后跳出该
结构；
若所有比较结果都为假，
则执行 otherwise 后面的一组
指令。
注意：与C语言不同，每个case语句后不要跟break语句。
clc
k = input('请输入一个小于10的正整数：');
switch k
case 1
disp('你输入的是符合要求的最小的数')
case {2,4,6,8}
用{ }构成的数组，称为 cell 数组。
disp('你输入的是符合要求的偶数')
case {3,5,7,9}
disp('你输入的是符合要求的奇数')
otherwise
disp('不符合输入要求！')
end
五、循环结构
1、 for循环
2、 while循环
3、 continue , break
1、 for 循环结构
for循环根据用户设定的条件，对结构中的命令反复执行
固定次数的操作，一般用于已知循环次数的情形。for循环的一
般格式如下：
for x = array
statement （循环体）
end
x 为循环变量，数组 array 的列数决定 for 循环的次
数。每次循环，x 依次取数组 array 的一列
在 for 后面的表达式中的数组 array 可以是任何合法
的MATLAB数组。
Matlab中i，j是虚数单位，若程序中涉及复数运算，一
定不能使用i，j作为循环变量。
for 循环举例 (1)
例：所谓水仙花数是一个3位数，其各位数字的立方和等于该数
本身。如：153 = 13 + 53 + 33 。请将所有的水仙花数输出至一
行数组，共有多少个水仙花数？
r=[];
for n = 100:999
a = fix(n/100);
b = fix(n/10) - a*10;
c = rem(n,10);
if n == a^3 + b^3 + c^3
r = [r,n];
end
end
disp(['水仙花数包括： ',num2str(r)]);
disp(['共有 ',num2str(numel(r)),' 个水仙花数']);
for 循环举例 (2)
使用for循环，绘图时，
先从第1点连线到第2点，然
后暂停，再从第2点连线到第
3点，再暂停，…，一直到最
后一点，这样可以实现最简
单的动画效果。
pause
pause(n)
pause on
pause off
clc
clear all
x = 0:0.1:2*pi;
y = sin(x);
axis([-1,7,-1.2,1.2]);
hold on;
for k =1:length(x)-1
plot(x([k,k+1]),y([k,k+1]));
pause(0.1);
end
使程序暂停，用户按下键盘上任意键后继续
使程序暂停n秒后，再继续执行
使后面的pause 语句起作用
使后面的pause语句不起作用
for 循环举例 (3)
Fibonacci Numbers:斐波那契数列
十三世纪的意大利数学家Fibonacci提出了一个有趣的问题：假定一对兔子出生满
两个月就可以生一对小兔子，之后每一个月又可以再生一对小兔子。假定现在有一对
刚生下来的小兔子，请问一年以后应该有几对兔子？
年初，只有1对小兔子。
一月，小兔子还没长大，所以还是只有1+0=1对。
二月，小兔子长成大兔子，开始生下一对小兔子，共有1+1＝2对。
三月，大兔子又生了一对小兔子，而小兔子还没长大，所以共有2+1=3对。
四月，第一对小兔子也长大开始生小兔子，这个月生了两对，总共有3+2=5对。
五月，第二对小兔子也长大了，所以有三对大兔子会生小兔子，总共有5+3＝8对
。
六月，第四月生的小兔子也长大了，所以这个月生了5对，总共有8+5＝13对。
… ...
f1  1,
f 2  2,
f n  f n 1  f n  2
 n  3
n = 12;
f = zeros(n,1);
f(1) = 1;
f(2) = 2;
for k = 3:n
f(k) = f(k-1) + f(k-2);
end
f
f=
1
2
3
5
8
13
21
34
55
89
144
233
2、 while 循环结构
while 循环一般用于不知道循环次数的情形。while 循
环的一般格式如下：
while expression （条件）
statement（循环体）
end
expression的值为逻辑真（非0），则执行循环
体，直到表达式值为假
一般情况下，表达式 expression 都是标量，但
MATLAB 允许它为数组，此时只有数组元素都为真时，循环
体才被执行。如果表达式为空数组，被认为是假
while 循环举例 (1)
reply = 'No'
while ~strcmp(reply,'YeS')
reply = input('输入YeS，不然很烦的。 ', 's');
end
v = input(‘str’,’s’) 用这种形式使用input()函数，用户的键盘输
入数据被认为是字符，创建的变量 v 为字符型数据。输入数字
1，认为是字符‘1’，输入字母a，认为是字符‘a’
strcmp(s1,s2)
比较两个字符串s1,s2是否相同，相同则
返回逻辑真(0)，否则返回逻辑假(0)
while 循环举例 (2)
例：一张纸厚 0.06 mm 且足够大，试问将纸对折多少次，其
厚度将超过10000 m？
h = 0.06e-3
n = 0;
while h < 1e4
n = n + 1;
h = 2*h;
end
disp(['对折 ',num2str(n),' 次其厚度将超过10000米'])
while 循环举例 (3)
The 3n+1 Sequence
a famous unsolved problem in number theory.
Start with any positive integer n. Repeat the following steps:
 If n = 1, stop.
 If n is even, replace it with n/2.
 If n is odd, replace it with 3n + 1.
For example, starting with n = 7 produces：
7, 22, 11, 34, 17, 52, 26, 13, 40, 20, 10, 5, 16, 8, 4, 2, 1.
n = 7;
y = n;
while n > 1
if rem(n,2)==0
n = n/2;
else
n = 3*n+1;
end
y = [y n];
end
disp(y)
3、 continue , break
continue 和 break 一般与 if 语句配合使用
continue ：退出本次循环，执行下一次循环
break ：终止当前的while、for循环，跳至相应的end后面的语句
下面这段代码的输出是什么？
for k = 1:3
for m = 1:3
if m == 2
continue
end
disp([k,m])
end
end
这段代码的输出又是什么？
for k = 1:3
for m = 1:3
if m == 2
break
end
disp([k,m])
end
end
统计magic.m文件中的代码行的行数,不计空行和注释行
（使用continue命令）。
fid = fopen('magic.m','r');
count = 0;
while ~feof(fid)
line = fgetl(fid);
if isempty(line) || strncmp(line,'%',1)
|| ~ischar(line)
continue
end
count = count + 1;
end
fprintf('%d lines\n',count);
fclose(fid);
可edit magic.m打开文件查看，若统计结果不正确，该
如何修改上面的代码？
使用while循环按行读取fft.m文件的内容至一字符数组。
当遇到第一个空行时，用break退出while循环。（第一个空
行前的内容为函数fft的使用帮助）
fid = fopen('fft.m','r');
s = '';
while ~feof(fid)
line = fgetl(fid);
if isempty(line) || ~ischar(line)
break
end
s = sprintf('%s%s\n', s, line);
end
disp(s);
fclose(fid);
4、程序终止语句：return
程序代码一般而按流程执行完毕正常退出，但当遇到
某些特殊情况，程序需要立即退出时，就可以用return提前终
止程序运行。return语句结束所在脚本(函数)的执行，返回
command window(调用函数)。
clc
clear all
for k=1:5
if k==3
continue
end
disp(k)
end
disp(k)
clc
clear all
for k=1:5
if k==3
break
end
disp(k)
end
disp(k)
clc
clear all
for k=1:5
if k==3
return
end
disp(k)
end
disp(k)
分析一下这三段程序的输出分别是什么？上机验证。
循环练习
1、找出所有这样的 4 位正整数，其各位数字的4次方的和
等于该数本身，例如：1634 = 14 + 64 + 34 + 44
输出所有的满足条件的数至一行数组，这样的数共有多少个？
2、猴子吃桃问题： 有一堆桃子不知数目，猴子第一天吃
掉一半，觉得不过瘾，又多吃了一个，第二天照此办理，吃掉
剩下桃子的一半另加一个，天天如此，到第十天早上，猴子发
现只剩一只桃子了，问这堆桃子原来有多少个？
3、假设matlab没有提供sum( )函数，也没有数组运算功能，
请重新编写一段matlab程序，使用循环，完成前面的线性拟合
的例子。
六、错误控制结构：try - catch - end
在程序设计中，有时候会遇到不能确定某段代码是否
会出现运行错误的情况。这时候可以使用错误控制结构。
程序运行时，首先尝试执行try和catch之间的代码段，
如果代码执行没有错误发生，则程序通过，不执行catch和
end之间的部分，而是继续执行end后面的代码。
一旦try和catch之间的代码执行发生错误，则立即转而
执行catch和end之间的部分，然后才继续执行end后的代码。
matlab提供了lasterr函数，可以获取出错的原因。
try
statement1
总被执行,正确则跳至end后
catch
statement2
statement1出错，才执行这里
end
例：引用数组x中的元素，当“下标”不合法时，改为引用数
组的最后一个元素。
clc
重复运行这段程序，
clear all
每次输入一个不同的数。若
x = -10:3:10
输入不是1到7的整数，catch
k = input('请输入一个数：');
后的语句就会执行。
try
disp(x(k));
disp(['上面显示的是数组x的第',num2str(k),'个元素。']);
catch
disp(x(end));
disp(['第',num2str(k),'个元素不存在。try内的语句出错啦，上
面显示的是数组x的最后一个元素。']);
end