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