Matlab速记

发表于

Matlab笔记

前言——最终成绩构成

平时30+期中20+期末50

↓

平时包括： 平时作业与单元检测20+云班测试5+mooc单元检测5

命令行的常用快捷指令

clear 清楚工作空间中的变量

clear v1 v2 （v1 v2是变量名) 清除指定变量

clc 清楚命令行窗口内容

help 命令名 给出指定命令用法介绍

doc 命令名 打开指定命令帮助文档

whos 列出当前所有变量详细信息

↑↓ 用于浏览历史记录 双击执行

常用赋值语句四格式

表达式

变量名 = 表达式

变量名 = 函数名（参数）

[变量名列表] = 函数名（参数）

矩阵相关及其用法

magic （幻方） M = magic(n)

空矩阵 [ ]产生一个空矩阵，或者删除元素用

A(i) 列拼接后第i个元素

A(:) 列拼接

A(a:b:c)

A(x,y)第x行第y列

tril(A)下三角

triu(A)上三角

diag(A)取对角

diag(diag(A))

解线性方程组Ax = b

A\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
>>> A= rand(3,3)

A =

    0.8147    0.9134    0.2785
    0.9058    0.6324    0.5469
    0.1270    0.0975    0.9575

>> x= rand(3,1)

x =

    0.9649
    0.1576
    0.9706

>> b = A*x

b =

    1.2004
    1.5045
    1.0673

>> x_compa = A\b

x_compa =

    0.9649
    0.1576
    0.9706

>>

inv(A)*b

A^(-1)*b

还有一种课外的方法

rref [A b]（线性代数最后一列是解向量）

1
2
3
4
5
6
7
8
9
>> rref([A b])

ans =

    1.0000         0         0    0.9649
         0    1.0000         0    0.1576
         0         0    1.0000    0.9706

>>

矩阵特征值

用法/函数

lambda = dig（A）求特征值

dig 既可以找出特征值，也可以找出特征向量，使得计算变得容易
A^n 在matlab已经优化

例子

减肥食谱中的线性代数
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
>> A = [36 51 13;52 34 74; 0 7 11]

A =

    36    51    13
    52    34    74
     0     7    11

>> b = [33;45;3]

b =

    33
    45
     3

x =

    0.5369
    0.2370
    0.1219
Leslie模型及其应用

简单迁移模型：每年甲镇人口10%前往乙镇，乙15%前往甲。假设某年甲乙两镇人口各有120 80人，问两年后两镇人口数量分布如何？

设总人口不变，人口只在甲乙间流动

引入变量x1（k），x2（k）

有如下关系：

1
2
3
4
5
6
7
8
x1（k+1） = 0.9x1(k)+0.15x2(k)
x2 (k+1)  = 0.1x2(k)+0.85x2(k)

A = [0.9 0.15;0.1 0.85]

X(k) = [x1;x2]

X(0) = [120;80]

可知有这样的关系：X(k+1) = AX(k)

两年后（二阶）即求A(2)

Q:那长期趋势？在无限长时间后两镇人数的变化规律与矩阵A有何关联？

编辑器的常用快捷指令

F5 保存并运行当前脚本

ctrl+r 对选定的行做注释

ctrl+t （对应）取消注释

ctrl+] 对选定的行缩进（可重复使用）

ctrl+[ 对选定的行取消缩进

创建数组

逗号操作符 等差数列 公差 项数

冒号操作符 等差数列 末项 公差

linspace(a,b,n)

X,Y = meshgrid(x,y);%1：生成平面网格坐标矩阵 meshgrid 可以快速生成坐标数对

拼接矩阵

横向拼接 空格 要求列相同

纵向拼接 分号 要求行相同

符号计算

1
2


sym x

syms x y

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

### simplify() 与 expand()

### 符号表达式中变换替换方法——subs()

```matlab
subs(s，old,new)%将s表达式中的old变量替换成new
subs(s)
subs(s,new)
%或者
syms x y z
f = x^2+y^2+z^2
fval = subs(f,[x,y],[1,2])%[old1,old2][new1 new2]或者
fval = subs(f,{'x','y'},{1,2})%或者
fval = subs(f,{x,y},{1,2})

old可以是单一变量，也可以是多变量构成的向量，new是用来替换的符号表达式

符号计算精度及其数据类型转换

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

```digits(n)```设置vpa计算结果的有效数字的位数

```vpa(s)```计算符号表达式s的数值结果

```double(s)```将符号表达式s转化为双精度数值/

```char(s)```将符号表达式s转化为字符串

```vpa(s,n)```采用n位有效数字计算精度求s的数值结果

eg：

```matlab
ans = pi;
vpa(ans,10);
___________________________________________
ans =

3.141592654

级数求和

1
S = symsum(f,n,a,b)%级数n，f是表达式，a是起始点，b是终点（无穷级数可以用inf）

taylor 泰勒多项式（包含麦克劳林）

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
taylor（f,n,x) %n-1 mkll

taylor(f,n,x,a) %point a mkll

%high version:

taylor(f,x)%默认'Order'=6，5阶麦克劳林多项式展开

taylor(f,x,a)%a点泰勒公式展开

taylor(f,x,'ExpansionPoint',a,'Order',n)%n-1阶泰勒展开在a点

taylor (f，[x,y,z])%多元泰勒展开

taylorf(,[x,y],[1,1],n)%n阶分别对x y （1，1）展开

牢记n-1

微积分符号计算

求导

1
2
3
4
5
diff(f)%（不推荐）

diff(f,v)

diff(f,v,n)

积分

1
2
3
4
5
int(f)%（不推荐）

int(f,v)

int(f,v,a,b)

intlin

作图

plot

1
```plot(x1,y1,s,x2,y2,s,...)

多条可以一条语句表达

ezplot（fun [min max min max]）填入表达式以及约束条件即可
eg：绘制衰减震荡曲线

fplot可替代前者（2016b版本已经移除前者）

1
2
3
y = @(x) exp(-0.5*x).*sin(5*x);
fplot(y)
title('衰减震荡函数')

cylinder

cylinder - 生成圆柱，此 MATLAB 函数 返回半径等于 1 的圆柱的 x、y 和 z 坐标。该圆柱绕其周长有 20 个等距点。

1
2
3
4
5
cylinder([1 1])%上圆下圆半径都是1 未加第三个参数则默认20个矩形组成侧面

cylinder([1 1],3)%第三个参数代表由三个矩形组成侧面

cylinder([0 1 0],50)%分成上下两部分对称锥体

图形注解

legend

title(‘’) 添加标题

xlabel 调整x轴属性

ylabel调整y轴属性

数组运算加“.”

end 的特殊用法

非线性方程求跟与曲线拟合

非线性方程：f（x) = 0 有代数方程与含三角函数，指数函数或超越函数两种

下看一些常用的矩阵运算函数命令

solove命令

试着求解一个方程（组）

1
2
3
4
5
6
7
8
9
10
11
>> syms x
>> eqn = (x+2)^x ==2;
>> s = solve(eqn,x)
警告: Cannot solve symbolically. Returning a numeric approximation instead.
> In solve (line 303) 
 
s =
 
0.69829942170241042826920133106081
 
>> `

roots命令

多项式方程求根

p是次数由高到低排列的的多项式系数

1
2
3
4
5
6
7
p = [3 -2 4]%代表3x^2-2x-4 = 0
r = roots(p)

r =

   0.3333 + 1.1055i%虚数单位，所以命名变量注意避免重复冲突
   0.3333 - 1.1055i

用数值方法求解非线性方程—————常先画图估计根大概所在区间，之后逼近收敛

先找出隔根区间，再从一点或多点出发，逐次逼近，寻求满足精度的根的近似值

求根：二分法

求根： 不动点迭代法

定理1：设x*为x = @（X）的不动点，@’（x）的某邻域连续，且绝对值小于一，这不动点迭代法xk+1 = @（xk）局部收敛

定理2：…

Newton迭代法

曲线拟合

应用：人口增长问题拟合

微分方程

求解微分方程的欧拉法

CLASSIC Example

输出斐波那契数列

循环ver

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
function [ F ] = feiboNum( n )
%默认输入的n大于等于三
F = [1 1];
for k = 3:n
    F(k) = F(k-1)+F(k-2);
end
end

———————————————————————————————
%命令行输出如下
>> feiboNum(10)

ans =

     1     1     2     3     5     8    13    21    34    55

递归ver

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

function [ F ] = feiboDigui( n )
% 应用递归不能n 默认考虑大于 2 的数
if n<= 2
    F = [1,1];
else
    F = feiboDigui(n-1);
    F = [F,F(n-1)+F(n-2)];
end
_______________________________
>> feiboDigui(10)

ans =

     1     1     2     3     5     8    13    21    34    55

输入一个任意位整数数，单独输出他每个位数整数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
function [ digits ] = getDigits( n )
%n = input ('n=')不允许，函数直接f（n）形式使用
leastDigit = mod(n,10); %rem(n,10)
digits = leastDigit;
n = (n - leastDigit)/10; %fix(n/10)
while n ~=0
   leastDigit = mod(n,10); %rem(n,10)
   digits = [leastDigit,digits];
   n = (n - leastDigit)/10; %fix(n/10) 
end
end

_______________________________
>> getDigits(1254)

ans =

     1     2     5     4

>>

输出x^-2y = 1所有正整数解x,y<=1000

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 [ M,L ] = findSolution% 无输入
y = (1:1000)'; %也可 3：43
x = sqrt(2*y + 1);
id = find(fix(x)==x);%找取整
M = [x(id),y(id)];
L = size (M ,1);%行数
end
_______________________________
>> findSolution

ans =

     3     4
     5    12
     7    24
     9    40
    11    60
    13    84
    15   112
    17   144
    19   180
    21   220
    23   264
    25   312
    27   364
    29   420
    31   480
    33   544
    35   612
    37   684
    39   760
    41   840
    43   924

ps: 改进可以更matlab 避免循环

输出0到任意正整数范围内的质数

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
function [ p ] = myPrimes( n )
% 2, 3, ......,n % 偶数不考虑
p = 2;
for  k = 3:2:n
    %如果 k 是素数
    if isPrime (k)
        p = [p,k];
    end
end
    function tf = isPrime( k )      %自动缩进，子函数
        for f = 3:2:sqrt(k)
            if mod(k,f) == 0
                tf = 0;
                return;
            end
        end
        tf = 1;
    end
end
_______________________________
>> myPrimes(40)

ans =

     2     3     5     7    11    13    17    19    23    29    31    37

比较不同算法实现同一需求的运行速度 ‘cputime’

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

clc
t0 = cputime;
feiboNum(50000);
ts1 = cputime - t0

t0 = cputime;
feiboDigui(50000);
ts2 = cputime - t0
_______________________________
ts1 =

     0


ts2 =

    1.1875

输入整数n，用函数返回n层杨辉三角

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
function Y=young(n)
%默认返回五阶
if nargin==0,n=5;end
Y=eye(n);%以单位矩阵开始计算
Y(：，1)=ones(n,1); %第一列全变成1
for k = 3:n %k是行指标
    Y(k,2:k) = Y(k-1,1:k-1)+Y(k-1,2:k);%eg 第三行为例，为前面11000中的110 加上011成为121
end
end
——————————————————————————————————————————
>> young(5)

ans =

     1     0     0     0     0
     1     1     0     0     0
     1     2     1     0     0
     1     3     3     1     0
     1     4     6     4     1

分段函数

1
2
3
4
5
6
7
function f = mPiecewise2(x)
f = (x<0).*(x.^2)...%利用布尔代数代替if else
    +(x>=0&x<10).*()...
    +().*()


end

或者

1
2
3
4
5
6
7
8
9
10
11
12
function f = myPieceWise(x)
L = length(x);%向量对应分量个数 当然可以同理用于矩阵 L = length（x(:)), 为拉伸操作
for k = 1:L
    f(k) = piece(x(k));%向量分成各个分量

end
function f = piece(x)
%支持标量计算
if x<0
    f = x.^2;
elseif x>=0&x<10
    f = ...

汉诺塔

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
function hanoi2(n,A,B,C)

%strcat 改成sprintf 更容易控制格式
if nargin==1, A='A';B='B';C='C';end
if n==1
    disp(sprintf('No %d: %s -> %s',n,A,C))
else
    hanoi2(n-1,A,C,B)
    disp(sprintf('No %d: %s -> %s',n,A,C))
    hanoi2(n-1,A,C,B)
end
____________________________________________

>> hanoi2(3, 'A', 'B', 'C')
No 1: A -> C
No 2: A -> B
No 1: A -> C
No 3: A -> C
No 1: A -> C
No 2: A -> B
No 1: A -> C
>> hanoi2(3, '甲', '乙', '丙')
No 1: 甲 -> 丙
No 2: 甲 -> 乙
No 1: 甲 -> 丙
No 3: 甲 -> 丙
No 1: 甲 -> 丙
No 2: 甲 -> 乙
No 1: 甲 -> 丙

作图 旋转的内接正方形

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
points = [1 -1 -1 1 1
          1 1 -1 -1 1];%创建正交矩阵
x=  points(1,:);
y = points(2,:);
plot (x, y),hold on
axis equal off
pause(0.1)
ang = 12/180*pi
r = 1/(cos(ang) + sin(ang));
rotM = r*[cos(ang),-sin(ang)
         sin(ang),cos(ang)];
for k = 1:20
    points = rotM * points;
    x = points(1,:);
    y = points(2,:);
    plot(x,y)
end

作图 n叶玫瑰线 (polar实例)

1
2
3
4
5
%三叶玫瑰线
%方程 ρ = acos3θ (n=3) n决定叶数 n越大 就现需要theta所需要周期越大
theta = 0:0.001:4*pi;
r = cos(2.5*theta);
polar(theta,r,'r')

作图 条形图/直方图绘制 （bar bar3实例）

1
2
3
4
5
6
7
8
9
10
%条形图/直方图的绘制
%bar (2D)  bar3 (3D)
x = -2.9:0.2:2.9;
y = exp(-x.*x);
subplot(1,2,1);
bar(x,y)
title('../figure:2-D Bar Chart')
subplot(1,2,2)
bar3(x,y)
title('../figure:3-D Bar Chart')

作图 二元函数绘图 （plot3实例）

1
2
3
4
5
6
7
8
%二元函数绘图
%空间曲线的绘制：plot3（x1,y1,z1,S1,x2,y2,z2,S2)
%三维螺旋线： x = sin(t),y = cos(t), z = t; t是渐增的
t = linspace (0,10*pi,1000)
plot3(sin(t),cos(t),t,'linewidth',2)
%xlabel('sin(t)','fontsize',18)
%ylabel
%zlabel 皆是参数调整 不作要求 可以生成后用属性编辑器鼠标可视化操作

作图 空间曲面绘制 生成地球经纬球像（mesh系列实例）

1
2
3
4
5
6
7
8
9
10
11
%空间曲面绘制的基本步骤 
%X = 1:6;
%Y = 1:8;
%X,Y = meshgrid(x,y);%1：生成平面网格坐标矩阵 meshgrid 可以快速生成坐标值
%Z = f(X,Y);%2: 计算网格点上的函数值 是数组运算决定
% x y z 都是 α β γ坐标角度
%mesh(X,Y,Z);%%只显示网格面
%surf(X,Y,Z);%%填充中间区域
%neshc()生成等高线
%meshz()成平台
%3 绘制网面
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
%下为mesh（）实例1
function earth
t = linspace(-pi/2,pi/2,20);%纬度数组
theta = linspace(0,2*pi,50);%经度数组

%上面已经生成20*50个数准备对应

[t1,t2] = meshgrid(t,theta);

%对应成1000个点

x = cos(t1).*cos(t2)%点乘是数组对应位置元素相乘，俩数组必同型
y = cos(t1).*sin(t2)
z = sin(t1);
mesh(x,y,z);
colormap('Cool');
axis square off%椭球修正成球
end

作图 等高线以及立体等高线绘制 （contour contour3实例）

1
2
3
4
5
6
7
%contour 与 contour3
subplot(1,2,1)
contour(x,y,z,20);%多了个参数 为z坐标的取值数目 即取值区间
title('')
subpolot
colorbar%含刻度的色相比照表
contour3(x,y,z,20);

作图 点线连接

X = [x1,x2,…]
Y同

1
2
3
4
5
6
7
X = [2,2,2,2];
Y = [0,2,2,0];
Z = [0,0,4,0];%点(2,0,0)(2,2,0)(2,2,4)组成三角形（最后闭合）
A = [1,0,0,1];
B = [0,1,0,0];
C = [0,1,2,0];%同理
plot3(A,B,C,X,Y,Z,'-');

创建符号数组f（2行2列）

使用syms 定义20个符号变量x1 to x20（for i = 1：20）

1
2
3
for i = 1:20
    eval(sprintf('syms x%d',i))
end

计算 f = arctan(x) 的导数并绘图

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
syms x
y = atan(x);
y1 = diff(y,x)
y2 = diff(y1,x)
figure,fplot(y)
figure,fplot(y1)
figure,fplot(y2)
____________________________________________
y1 =

1/(x^2 + 1)

y2 =

-(2*x)/(x^2 + 1)^2

符号表达式中变换替换方法——subs()对应例子

1
2
3
4
5
6
7
8
9
10
11
t = linspace(-6,6,30);%-6到6等距离取了30个点
r = 3-t.^2/36;
cylinder(r,30);
axis off%隐藏坐标系
colormap('jet')
syms y h
f = 3 - y*y/36;%被积函数
V = pi* int(f*f,y,-6,h)%积分表达式，依次是目标函数，积分变量，积分下限，积分上限
V12 = subs(V,h,6)%替代 把V中的h用6代替 ans放入V12中 用来计算某些积分表达式

%%%12米深啤酒桶形状，顶部底部半径2米，中心半径3米。

数值转化为数值数据：double(A)实例

练习：计算曲线段f(X) = exp(ax)sin(bx),0=<x=<2pi,绕x轴旋转曲面体积

1
2
3
4
5
6
syms a b x
f = exp(a*x)*sin(b*x);
f1 = subs(f,a,-0.2);
f2 = subs(f1,b,0.5);
V = pi*int(f2*g2,x,0,2*pi)
double(V)

泰勒实例：椭圆积分近似表达式

$\int_{0}^{\frac{\pi}{2}}$ $\sqrt{1-e^2*cos^2(t)}$dt 

1
2
3
4
5
6
syms t e2 x
f = sqrt(1-e2*x);
F = taylor(f,x,0,'Order',6)
g = gubs(F,x,cost(t)^2);
s = int(g,t,0,p/2);
double(subs(s,e2,0.8))

微积分实验

一

已知函数y = a*exp(x)/sqrt(a^2+x^2) 求x = 5a处二阶导数值

二

六阶mkll多项式 y= 1/1+x^2

求定积分 $\int_{？}^{\frac{？}{？}}$ $\sqrt{1-e^2*cos^2(t)}$dt

…

四

funy = 174.42(x1/x5)(x3/(x2 - x1))^0.85*
sqrt((1 - 2.62(1 - 0.36(x4/x2)^(-0.56))^(3/2)(x4/x2)^1.16)/(x6x7));

周一课上实验：线性代数实验

动物年龄组的变化趋势

三个年龄组：零到五岁，六到十岁，十一到十五岁。

第二组平均繁殖四个后代

第三组平均繁殖三个后代

第一第二组五年的存活了吧分别为0.5，0.25；

现在又三个年龄组动物各1000；计算五年后，十年后，十五年后各个年龄组动物数量

模型构建

t0=0；t1=5,t2=10,t3-15

各个年龄组动物数量关系