Lec 1：MATLAB 入门、数组、数值精度

1. 课程基本信息

1.1 课程名称

Introduction to Scientific Computing 科学计算导论

本课程围绕 MATLAB 展开，学习科学计算中常用的编程、矩阵运算、可视化、数值方法和工程应用。

2. Course outline & schedule 课程安排

3. 参考资料

3.1 Textbook 教材

Essential MATLAB for Engineers and Scientists 作者：Hahn & Valentine

3.2 其他资料

• LLM：GPT、Gemini、DeepSeek 等
• MATLAB Help：MATLAB help 可能是最好的教材
• Course website：https://course.zju.edu.cn/

4. 课程内容简介

4.1 What is this course about?

本课程涉及三个方面：

1. About computer programming

编程语言就是给计算机的指令。

MATLAB 是：

• 一种高级编程语言
• 一个专门用于工程和科学计算的计算系统
• 特别适合：
  • 数值计算
  • 数据分析
  • 可视化

2. About scientific computing

科学计算关注：

• 解决科学问题的基本技术
• 解决科学问题的策略
• 应用领域：
  • 科学
  • 工程
  • 经济学
  • 其他很多学科

3. About programming methodology

编程方法论包括：

• 程序设计
• 编程风格
• 调试 debugging
• 算法和数据结构
• 抽象的计算机科学 / 编程概念

5. Course topics 课程主题

本课程会涉及：

• MATLAB 基础编程
• MATLAB 中的矩阵 / 数组操作
• 可视化 visualization
• 线性系统 linear systems
• 调试代码 debugging
• 模拟动态系统 simulating dynamic systems
• 优化问题 optimization
• 组合问题 combinatoric problems
• 随机化 randomization
• 数据分析 data analysis
• 信号处理 signal processing
• 控制 control

6. Course mechanics 课程机制

6.1 上课时间

• Lectures：8 周，每周 3 小时
• Tuesday：13:25 – 15:50
• Room：西 2-310

6.2 成绩组成

Final project 要展示你的 MATLAB 能力。

6.3 作业提交

• 提交时间：Wednesday 1:00 pm，课前
• 迟交政策：
  • 迟交 48 小时内：扣 25%
  • 迟交 72 小时内：扣 50%
  • 超过 72 小时：不接受

6.4 Plagiarism policy 抄袭政策

• 作业 / 项目可个人完成或 pair programming。
• 如果发现作业整体或部分复制，所有相关人员该次作业为 0。
• 第二次违规：
  • 课程直接不及格
  • 上报相关部门
• 代码会进行结构相似性分析。
• 仅修改变量名、加几行代码，不能规避抄袭检测。

7. MATLAB 简介

7.1 MATLAB 是什么？

MATLAB = MATrix LABoratory

课程主要使用 MATLAB。

PPT 中提到的版本：

• R2019b
• R2013a
• R2012a
• 等

浙江大学相关访问地址：

1
http://ms.zju.edu.cn/

8. 为什么学习科学计算？

科学计算可用于：

8.1 Simulation 模拟

模拟物理、化学、生物、环境、社会系统。

现实实验可能：

• 昂贵
• 缓慢
• 危险
• 不可能实施

因此可以用计算机模拟部分替代真实实验。

8.2 Visualization 可视化

用计算机可视化复杂现象，例如：

• Rayleigh-Taylor instability
• 两种流体混合
• 复杂动态系统

8.3 Analysis 分析

计算系统性质：

• 应力 stresses
• 温度 temperatures
• 质量分布 mass distributions
• 浓度 concentrations

8.4 Optimization 设计优化

通过计算找最优设计参数。

8.5 Nanoscale and Atomistic Systems

分析纳米尺度原子相互作用，从而推导材料性质。

8.6 Combinatorics 组合问题

在大量可能性中寻找答案。

8.7 Data Analysis 数据分析

科学和商业越来越依赖数据。 通过分析大量数据：

• 发现趋势
• 发现模式
• 验证假设

例如 Human Genome Project 人类基因组计划。

8.8 Systems Control 系统控制

现代系统设计中，计算机常作为中央控制单元：

• 从传感器不断测量数据
• 设置控制输出

9. Scientific computation 科学计算

科学家和工程师会：

1. 建立物理问题的数学模型
2. 计算数学模型的解
3. 求解非线性方程
4. 求解优化问题
5. 求解微分方程，以理解真实世界如何动态变化

9.1 为什么需要 simulations？

真实实验可能：

• expensive 昂贵
• slow 慢
• perilous 危险

计算机模拟可替代或部分替代真实实验。

10. Approximation 近似

科学和工程模型通常都包含近似。

对于物理系统：

• 所有模型都包含近似
• 有些模型在一定容忍误差内仍然有预测能力

例子：

• Newton’s laws 牛顿定律
• Relativity 相对论

10.1 近似相关问题

需要考虑：

• Model errors 模型误差
• Solution errors 解的误差
• 近似的适用范围
• 需要的精度
• 如何验证计算结果 verify calculations
• 如何验证模型 validate model

11. Aspects of computation 计算的两个方面

11.1 Hardware 硬件

运行计算的机器。

包括：

• RAM
• CPU
• power supply
• hard drive
• motherboard
• videocard

11.2 Software 软件

指导计算的程序。

11.3 计算机例子

ENIAC

• 日期：1947/7/29
• 第一台计算机
• 5000 flops
• 用于弹道轨迹计算

Blue Gene

• 时间：2012/06
• Total Cores：1,572,864
• Power：7890 KW
• Memory：1,572,864 GB
• Rmax：16,324.8 Tflops
• Rpeak：20,132.7 Tflops

Human brain

估计约比 Blue Gene 强 5 倍。

12. Machine Model 机器模型

12.1 计算机如何处理信息？

计算机系统包括：

• 高密度存储：hard drive
• 随机访问内存：RAM
• 主板和 CPU
• 外设 / 用户界面
• 电源

12.2 从高级代码到硬件

程序执行的大致层次：

1
Applications
2
High level code
3
Assembly / machine code
4
Voltages
5
Microprocessor
6
Individual transistors

也就是说，我们写的 MATLAB 代码最终都会被转换成底层机器能够执行的形式。

13. Assembly language 汇编语言示例

如果直接在机器 / 汇编层编程，会很麻烦。

13.1 伪指令集

13.2 高级语言的意义

汇编语言编程很痛苦，因此出现了：

• Java
• C
• MATLAB

以及：

• interpreters 解释器
• compilers 编译器

它们把人类可读代码转换成低级指令。

例如高级代码：

1
sum = 0;
2
for i = 1:500
3
    sum = sum + x(i);
4
end
5
avgVal = sum/500;

会被编译 / 解释成底层机器代码。

14. 为什么是 MATLAB？

MATLAB 部分源自 FORTRAN。

FORTRAN = FORmula TRANslation 最初就是为数学公式计算设计的。

MATLAB 使复杂数学表达式求值非常方便。

15. MATLAB Help 和 File Exchange

15.1 MATLAB help system

需要熟悉 MATLAB help 系统。 它包含所有 MATLAB 函数的文档。

常用命令：

1
help 函数名
2
doc 函数名
3
lookfor 关键词

15.2 MATLAB File Exchange

提供大量预写代码和示例。

16. MATLAB 界面

MATLAB 界面通常包括：

17. MATLAB Statements 语句

MATLAB 最基本的命令是 statement。

例：

1
foo = x^2 + sin(5*y) / exp(67*z);

含义：

• 右边是表达式 expression
• 左边是变量名
• MATLAB 计算右边结果并赋值给左边变量

17.1 分号

1
a = 5;

分号的作用：抑制输出。 如果没有分号：

1
a = 5

MATLAB 会在命令窗口显示：

1
a =
2
     5

18. MATLAB Built-in Functions 内置函数

MATLAB 有大量内置函数，例如：

1
sin(x)
2
cos(x)
3
exp(x)
4
log(x)
5
atan(x)
6
cosh(x)
7
sinh(x)
8
mean(x)
9
median(x)

几乎你能想到的数学 / 统计函数，MATLAB 都有。

19. Array basics 数组基础

MATLAB 中数组非常重要。 数组可以理解为一串数字。

19.1 创建数组

1
x = [1, 3, 9, 11, -10.2];

逗号 , 或空格都可以分隔同行元素：

1
x = [1 3 9 11 -10.2];

19.2 访问数组元素

MATLAB 下标从 1 开始，不是 0。

1
foo = x(3);

取出第三个元素。

19.3 修改数组元素

1
x(2) = 27;

把第二个元素改成 27。

20. Colon operator 冒号创建数组

20.1 基本格式

1
x = start:stop;

默认步长为 1。

例：

1
x = 1:15;

生成：

1
[1 2 3 ... 15]

20.2 指定步长

1
x = start:step:stop;

例：

1
x = 4:2:28;

生成：

1
[4 6 8 ... 28]

例：

1
x = -1.3:0.1:1.3;

从 -1.3 开始，每次加 0.1，直到不超过 1.3。

20.3 上界不会被超过

例：

1
x = 1:2:6;

结果：

1
x = [1 3 5]

不会包含 7，因为超过 6。

21. Operations on arrays 数组运算

21.1 对数组逐元素计算

1
x = 0:0.1:10;
2
y = exp(5*x).*sin(x);
3
plot(x,y);

注意 .*表示逐元素乘法。

21.2 数组加法

两个数组大小相同才能加。

1
x = [1 2 3 4 5];
2
y = [7 8 9 10 11];
3
z = x + y;

结果：

1
z = [8 10 12 14 16]

减法同理。

21.3 数组乘法和除法

逐元素乘法：

1
z = x.*y;

逐元素除法：

1
z = x./y;

补充： * 是矩阵乘法。 .* 是数组对应元素相乘。 / 是矩阵右除。 ./ 是对应元素相除。

21.4 内置函数通常支持数组输入

1
x = 1:10:1000;
2
y = sqrt(x);

sqrt(x) 会对 x 中每个元素开平方。

22. Scripts 脚本

22.1 什么是 script？

Script 是 .m 文件，包含一组 MATLAB 语句 / 命令。

22.2 创建脚本

命令行：

edit filename.m

22.3 运行脚本

如果文件名是 filename.m，运行：

filename

也可以按 F5 或点击 Run。

22.4 脚本作用域

脚本可以从：

• 当前目录 current directory
• MATLAB search path

中被找到并运行。

23. Current working directory 当前工作目录

MATLAB 有一个当前工作目录。

建议：

1. 创建一个专门放课程代码的文件夹。
2. 在 MATLAB 中把当前工作目录设置到那里。
3. MATLAB 会优先在当前目录查找脚本和函数。

24. Linear algebra primer 线性代数入门

线性方程组可以写成：

1
A * X = B

其中：

• A 是系数矩阵
• X 是未知向量
• B 是右端项

理论上：

1
X = A^(-1) B

MATLAB 中一般用：

1
X = A\B;

比显式 inv(A)*B 更推荐。

25. 打开、运行、中断 M-file

25.1 新建 / 打开脚本

edit newScript

25.2 执行脚本

• F5
• 命令窗口输入脚本名
• 工具栏 Run

25.3 中断程序

如果程序死循环或运行太久：

Ctrl + C

26. Memory 内存与二进制

26.1 bit 和 byte

• 计算机内存中每个位置存储有限个二进制位 bit。
• 常见存储大小：
  • 8 bits
  • 16 bits
  • 32 bits
  • 64 bits
• 8 bits = 1 byte

26.2 二进制表示非负整数

例：

1
1011₂ = 8 + 2 + 1 = 11₁₀
2
110010₂ = 32 + 16 + 2 = 50₁₀
3
11111111₂ = 255₁₀

8 bit 无符号整数最大是：

1
255

27. Overflow 溢出

内存位置大小固定，因此能表示的数值范围有限。

例：

1
X = uint8(78);
2
Y = uint8(190);
3
Z = X + Y;

理论上：

1
78 + 190 = 268

但 uint8 最大只能表示 255，因此结果会被截断 / 饱和到 255。

27.1 常见无符号整数范围

28. Floating point 浮点数

科学计算中不仅有整数，还有：

• 有理数，例如 0.335，-27.890
• 无理数，例如 pi 和 e

计算机通常用有限精度近似表示它们。

28.1 科学计数法

例：

1
+2.1345e17

包含：

• sign 符号
• mantissa 尾数，例如 2.1345
• exponent 指数，例如 17

28.2 IEEE floating point

Single precision 单精度

总共 32 bits：

• 1 bit sign
• 8 bits exponent
• 23 bits mantissa

Double precision 双精度

总共 64 bits：

• 1 bit sign
• 11 bits exponent
• 52 bits mantissa

MATLAB 默认使用 double。

29. Floating point arithmetic 的问题

浮点数位数有限，因此每次数值运算都可能损失信息。

例如如果只能保留小数点后 2 位：

1
5.43 x 10^6
2
+ 7.43 x 10^4
3
= 5.5043 x 10^6

但存储时可能只能保存部分数字，后面数字会被舍入。

29.1 Take-home message

计算机并不像你想象的那样擅长精确算术。

浮点计算结果通常应看作近似值。

30. Numerical analysis 数值分析

为了获得可靠结果，需要仔细设计代码，避免或减少：

• roundoff error 舍入误差
• other imprecision 其他不精确性

专门研究这些问题的人叫：

1
numerical analysts

31. MATLAB number types

可以用 whos 查看变量存储类型。

常见类型：

| 类型 | 含义 | | — | | | uint8 | 无符号 8 bit 整数 | | int64 | 有符号 64 bit 整数 | | double | 双精度浮点数 | | single | 单精度浮点数 |

指定类型例子：

1
x = int32(67);

32. 数组与数值类型

数组中所有数值通常共享同一种数值类型。

1
x = uint16([1 3 4 5 89]);
2
x = single(0:0.1:100);

33. It’s all bits

计算机中所有信息最终都是 bits：

• 数字
• 文本
• 图像
• 声音

34. 文本、声音、图像表示

34.1 Text 文本

ASCII code 给每个字符和符号一个数字。

MATLAB 中字符串可看作整数数组。

MATLAB 用 16 bit Unicode 表示字符。

34.2 Sound 声音

声音可表示为随时间变化的采样数值数组。

34.3 Images 图像

图像是二维数值数组。

彩色图像每个像素存三个数：

• Red
• Green
• Blue

即 RGB 三个通道。

35. Useful MATLAB commands

例：

1
clear
2
clc
3
whos

Lec 2：Scripts、Functions、循环、条件、逻辑向量

1. M-file

M-file 是存储在 MATLAB .m 文件中的一组指令。

M-file 可以是：

• script 脚本
• function 函数

2. Functions vs. Scripts

2.1 Function 示例

1
function d = myFunFunction(a,b,c)
2
% Takes scalar values a, b, and c, modifies and sums the values
3

4
% Modifying the input values
5
var1 = 2*a;
6
var2 = b^2;
7
var3 = sqrt(c);
8

9
% Calculating the output
10
d = var1 + var2 + var3;

特点：

• 第一行以 function 开头
• 输入：a,b,c
• 输出：d
• 函数内部变量不会自动出现在外部 workspace

2.2 Script 示例

1
% Define values of scalars a, b, and c
2
a = 1;
3
b = 2;
4
c = 15;
5

6
% Modifying the scalar values a, b, and c
7
var1 = 2*a;
8
var2 = b^2;
9
var3 = sqrt(c);
10

11
% Calculating the sum of the modified scalar values
12
d = var1 + var2 + var3;

特点：

• 没有 function 开头
• 变量直接存在于当前 workspace 中
• 运行后 a,b,c,var1,var2,var3,d 都可能留在 workspace

3. Functions 和 Scripts 的共同点

二者都：

• 存在 .m 文件中
• 使用 Editor 窗口编辑
• 是独立代码块
• 可以从 Command Window 调用
• 可以被其他函数调用

4. 为什么用 Scripts？

4.1 Script 是什么？

MATLAB M-file，其中所有变量、常量等都在文件中调用。

4.2 为什么写 script？

因为：

• 可能需要多行代码解决问题
• 一直在 Command Window 重复输入很麻烦
• Command Window 不方便保存代码
• Script 便于注释、修改、debug

5. 为什么用 Functions？

5.1 Function 是什么？

以 function 开头的 MATLAB M-file。

通常：

• 传入 input
• 返回 one or several outputs

5.2 为什么使用 function？

函数是一个 black box 黑箱。

优点：

• 功能测试 / 验证后可反复使用
• 不必每次质疑输出
• 可把复杂计算拆分成小代码块
• 支持 top-down program design 自顶向下设计

6. Scripts vs. Functions 的形象理解

PPT 用“火箭”和“星球”类比：

1
function [out1,out2,out3] = mars(in1,in2)
2
function [out1,out2] = saturn(in1)

• 火箭携带 inputs 和 outputs。
• Mars 函数看不到 Earth 上的变量，除非变量作为 input 被火箭带过去。
• 函数内部有自己的 workspace。

7. Function 语法

7.1 单输入单输出

1
function output = myFunc(input)

7.2 多输入多输出

1
function [output1,output2,output3] = myFunc(input1,input2,input3)

注意：函数定义第一行应写在同一行。

7.3 命名和保存规则

• 函数名必须以字母开头。
• 保存文件名必须和函数名一致。

例如：

1
function y = squareMe(x)
2
    y = x^2;
3
end

文件必须保存为：

1
squareMe.m

8. Newton’s Method 牛顿法求平方根

目标：只用 + - * / 求给定数 a 的平方根。

考虑函数：

1
y(x) = x^2 - a

平方根对应：

1
x^2 - a = 0

即求这个函数的根。

8.1 牛顿法思想

在当前点 x：

1. 用切线近似函数
2. 求切线与 x 轴交点
3. 用交点作为新的 x
4. 重复直到足够接近

8.2 算法

初值：

1
x = a/2;

循环条件：

1
abs(x^2 - a) > 1.0e-6

每步：

1
y = x^2 - a;
2
gradient = 2*x;
3
step = -y/gradient;
4
x = x + step;

对应文件：

1
my_sqrt.m

函数：

1
my_sqrt(a)

9. 如何调用函数？

可以：

1. 从 command line 调用，必须包含 inputs。
2. 从其他函数调用。

例：

1
result = my_sqrt(2);

10. Function H1 line 和 help

10.1 H1 line

函数第一行注释称为 H1 line，可被 lookfor 搜索。

例：

1
function out1 = testFunction(in1)
2
% A test function
3
% This would report, in the help function, all of the functionality of the
4
% function
5
% in1 = function input (units: not specified)
6
% out2 = function output (units: not specified)

10.2 help 查询

函数开头连续注释块会作为 help 信息返回。

1
help testFunction

11. Recursive functions 递归函数

函数可以调用自己，这叫 recursive function。

例：

1
function recursiveAdd(nToAdd)
2
% Adds a number to the input up to
3
% a value of 25
4

5
if nToAdd < 25
6
    nToAdd = nToAdd + 1;
7
    disp(nToAdd)
8
    recursiveAdd(nToAdd);
9
end

递归必须有终止条件，否则会无限调用。

12. 函数调用时发生什么？

1. 实参 arguments 被复制到输入参数 input parameters。
2. 函数体在自己的 private workspace 中运行。
3. 函数不能看见调用环境中的变量。
4. 函数结束时，输出参数被复制回调用环境。
5. MATLAB 通过 search path 查找函数：
   • 当前目录优先
   • 然后搜索路径中的其他目录

13. Search paths 搜索路径

MATLAB 调用函数 / 脚本时，会搜索一系列文件夹，即 search path。

包含：

• current working folder
• MATLAB 安装时自带 toolbox 文件夹
• 用户添加的路径

13.1 查看路径

1
path

13.2 设置启动工作目录

1
userpath('C:\MyMATLABwork')

13.3 添加路径

1
addpath('某个文件夹路径')

13.4 删除路径

1
rmpath('某个文件夹路径')

14. 函数变量去哪了？

14.1 Script 变量会保留

除非清除，script 中产生的变量会留在 workspace。

14.2 Function 变量不会保留

函数内部变量局限于函数 workspace。

PPT 总结：

1
Whatever happens in a function, stays in a function.

15. 清理 workspace 和 command line

15.1 清除变量

1
clear

也可以清除指定变量：

1
clear x y

15.2 清空命令窗口

1
clc

15.3 注意

不要在 function 内随便调用 clear。

clear 通常在：

• script
• command line

中使用。

16. 注释代码

16.1 注释原则

一般规则：至少每三行代码有一行注释。

注释有助于：

• 别人理解代码
• 别人 debug
• 助教 / 老师评分

16.2 快捷键

17. Example 1：Monthly Loan Payment

固定利率按揭贷款。

输出：

• monthly payment c

输入：

• loan amount P0
• interest rate r
• number of monthly payments N

文件：

1
loan.m

18. Example 2：Evaluating an Expression

要求写一个函数文件计算某表达式，并允许 x 是向量。

计算：

• f(x) for x = 6
• f(x) for x = 1,3,5,7,9,11

文件：

1
Ex10_2.m

补充：如果函数要支持向量输入，应使用逐元素运算： .*、./、.^

19. Example 3：Converting Temperatures

写函数把华氏度 °F 转换为摄氏度 °C。

文件：

1
FtoC.m

公式：

1
C = 5*(F - 32)/9;

如果支持向量：

1
C = 5*(F - 32)./9;

19.1 热膨胀问题

长度变化：

1
DL = a*L*DT;

其中：

• a：热膨胀系数
• L：长度
• DT：温度变化

题目：

矩形铝板：

• 4.5 m × 2.25 m
• a = 23x10^-6 / °C
• 温度从 40°F 到 92°F

求面积变化。

20. Anonymous Function 匿名函数

适合短的一行计算。

20.1 一般形式

1
function_name = @(arguments) expression

例：

1
FtoC = @(F) 5*(F-32)./9;
2
cube = @(x) x^3;
3
circle = @(x,y) 16*x^2 + 9*y^2;

补充： 如果输入可能是向量，cube 应写成：

1
cube = @(x) x.^3;

20.2 匿名函数例子

1
FA = @(x) exp(x^2)/sqrt(x^2+5);
2
FA(2)

支持向量时：

1
FA = @(x) exp(x.^2)./sqrt(x.^2+5);
2
FA([1 0.5 2])

双变量：

1
HA = @(x,y) 2*x^2 - 4*x*y + y^2;
2
HA(2,3)

如果支持数组输入：

1
HA = @(x,y) 2*x.^2 - 4*x.*y + y.^2;

21. Subfunctions 子函数

一个 function file 可以包含多个用户定义函数。

• 第一个函数叫 primary function。
• 文件名对应 primary function。
• 其他函数叫 subfunctions。
• subfunctions 只在该文件内部可见。
• 每个函数都有自己的 workspace。
• 常用于 utility calculations。

22. Example 4：Mean and Standard Deviation

写函数计算一组数的：

• mean 平均值
• standard deviation 标准差

文件：

1
stat.m

23. Function Functions

MATLAB 内置函数 fzero 可找函数零点：

1
f(x) = 0

问题：如何把函数 f(x) 描述给 fzero？

两种方式：

1. Function handle
2. Function name string

24. Function Handles 函数句柄

函数句柄是一种数据类型，保存和函数相关的唯一值。

24.1 获得方式

1
@function_name
2
@cos
3
@FtoC

如果 FtoC 是匿名函数，则直接用变量名 FtoC。

例：

1
f = @cos;
2
f(0)

文件：

1
funplot.m
2
Fdemo.m

25. Function Names in Strings

较老方法，效率较低。

传入函数名字符串：

1
'cos'
2
'FtoC'

用 feval 执行：

1
var = feval('function_name', arguments);

文件：

1
funplotS.m

26. 解决真实问题的策略

1
Start simple and build to complexity

先从简单问题开始，再逐步增加复杂性。

27. 初始化向量

27.1 离散输入

1
x = [0.5, 1.5, 2, 3.2, 1.5, 0.4];

27.2 冒号

1
x = 1:1:5;
2
x = 1:5;
3
x = -5:2:5;

27.3 linspace 和 logspace

1
x = linspace(1,5,5);
2
x = logspace(1,5,5);

linspace(a,b,n)：从 a 到 b 均匀取 n 个点。 logspace(a,b,n)：从 10^a 到 10^b 按对数均匀取 n 个点。

28. 访问向量元素

假设：

1
x = [5, 9, 4, 1, 7, 3, 4, 8];

28.1 第一个元素

1
y = x(1);

结果：

1
y = 5

28.2 范围索引

1
y = x(1:3);

结果：

1
[5 9 4]

1
y = x(6:8);

结果：

1
[3 4 8]

28.3 最后一个元素

1
y = x(end);

结果：

1
8

28.4 偏移索引

1
y = x(1+4);

即 x(5)，结果：

1
7

28.5 倒数第三个元素

1
y = x(end-2);

结果：

1
3

28.6 离散索引

1
y = x([1 3 6]);

结果：

1
[5 4 3]

29. 修改 / 添加向量元素

假设：

1
x = [1, 2, 3];

29.1 覆盖第一个元素

1
x(1) = 5;

结果：

1
[5 2 3]

29.2 覆盖最后一个元素

1
x(end) = 1;

结果：

1
[5 2 1]

29.3 添加一个末尾元素

1
x(end+1) = 8;

结果：

1
[5 2 1 8]

29.4 添加多个末尾元素

错误写法：

1
x(end+1) = [6, 9];

因为左边只有一个位置，右边有两个数，维度不匹配。

正确写法：

1
x(end+1:end+2) = [6, 9];

结果：

1
[5 2 1 8 6 9]

30. 删除向量元素

用空数组 [] 删除。

假设：

1
x = [8, 4, 5, 9, 3, 2, 5, 6];

30.1 删除第一个元素

1
x(1) = [];

结果：

1
[4 5 9 3 2 5 6]

30.2 删除多个元素

1
x(end-3:end) = [];

结果：

1
[8 4 5 9]

31. for loops

公式向量化通常更快，但有时需要循环。

31.1 基本 for 循环

1
for jj = 1:20
2
    disp(jj)
3
end

输出 1 到 20。

当迭代次数已知时，for 循环优于 while 循环。

31.2 任意冒号表达式

1
for jj = 1:3:19
2
    disp(jj)
3
end

输出：

1
1 4 7 10 13 16 19

31.3 离散向量循环

1
for jj = [1 9 5]
2
    disp(jj)
3
end

依次输出 1、9、5。

31.4 嵌套循环

1
for jj = 1:10
2
    for kk = 1:10
3
        disp([jj kk])
4
    end
5
end

每一个 for 都需要一个 end。

32. Relational operators 关系运算符

关系运算符比较数值并返回 true 或 false。

MATLAB 中：

• true 表示 1
• false 表示 0

注意：

1
=

是赋值。

1
==

才是判断相等。

33. while loops

33.1 基本语法

1
while <relationalExpression>
2
    <expressionToEvaluate>
3
end

PPT 中 End 应为 MATLAB 语法 end，大小写一般 MATLAB 可接受，但习惯小写。

33.2 简单例子

1
a = 1;
2
while a < 10
3
    a = a + 1;
4
end

循环结束时：

1
a = 10

while 循环适合不知道迭代次数的情况。

34. for 和 while 实现 n 次迭代

34.1 for loop 牛顿法

1
n = 100;
2
a = 2;
3
x = a/2;
4

5
for ii = 1:n
6
    x = (x + a/x)/2;
7
end

34.2 while loop 牛顿法

1
n = 100;
2
a = 2;
3
x = a/2;
4
ctr = 0;
5

6
while ctr < n
7
    ctr = ctr + 1;
8
    x = (x + a/x)/2;
9
end

二者没有明显时间优势，取决于代码可读性和你习惯哪种。

35. Vectorization 向量化

考虑：

1
a = rand(1,N);

希望计算：

1
b_i = a_i * a_i

有三种方式：

1
a.^2
2
a.*a
3
for loop

向量化更快，通常可快几个数量级。

35.1 性能比较代码

1
clear; clc; close all
2

3
n = round(logspace(0,8,20));
4

5
for jj = 1:length(n)
6

7
    a = rand(1,n(jj));
8
    clear b
9
    tic
10
    for kk = 1:length(a)
11
        b(kk) = a(kk)*a(kk);
12
    end
13
    t_for(jj) = toc;
14

15
    clear b
16
    tic
17
    b = a.*a;
18
    t_vect(jj) = toc;
19

20
end
21

22
[ph] = plot(n,t_for,'k.-',n,t_vect,'r.-');
23
set(ph(1),'LineWidth',4,'MarkerSize',25)
24
set(ph(2),'LineWidth',4,'MarkerSize',25)
25
set(gca,'YTick',logspace(-10,10,21))
26
set(gca,'XTick',logspace(0,20,21))
27
set(gca,'XScale','log','YScale','log')
28
set(gca,'FontSize',20)
29
set(gca,'FontName','Arial')
30
lg = legend('For loop','Vector');
31
set(lg,'Location','NorthWest')
32
xlabel('Number of array elements')
33
ylabel('Evaluation Time (s)')

补充：

• tic 开始计时。
• toc 返回从上次 tic 开始经过的时间。

36. if statement

36.1 基本语法

1
if condition
2
    statements
3
end

例：

1
a = 1;
2

3
if a == 1
4
    disp('a is equal to 1! Yay!')
5
end

37. Nested if 嵌套 if

1
a = 1;
2
b = 2;
3

4
if a == 1
5
    if b == 2
6
        disp(['a is equal to 1',...
7
              'and b is equal to 2.',...
8
              'Double yay!'])
9
    end
10
end

38. Boolean expressions / logical operators

逻辑运算符可以避免过多嵌套。

38.1 用 && 合并条件

1
a = 1;
2
b = 2;
3

4
if a == 1 && b == 2
5
    disp(['a is equal to 1',...
6
          'and b is equal to 2.',...
7
          'Double yay!'])
8
end

39. if / elseif / else

当有多个条件时：

1
if condition1
2
    statements1
3
elseif condition2
4
    statements2
5
else
6
    statements3
7
end

PPT 例子提到 FEM shape function。

40. Logical operators 完整表

40.1 or 与 xor 真值表

区别：

• or：只要至少一个真就真。
• xor：只有一个真时才真，两个都真反而是假。

41. Logical NOT

NOT 只需要一个输入，返回逻辑相反值。

1
not(0)   % 1
2
~0       % 1
3
not(1)   % 0
4
not(3)   % 0

非零值在逻辑中都视为 true。

42. Logical scalars 逻辑标量

例：

1
a = 3 > 2;    % 1
2
a = 2 > 3;    % 0
3
a = 2 >= 2;   % 1
4
a = 2 ~= 2;   % 0
5
a = 2 == 2;   % 1
6
a = ~1;       % 0

逻辑类型比 double 更省内存：

• logical：1 byte
• double：8 bytes

43. Logical vectors 逻辑向量

逻辑向量是由 0 和 1 组成的向量，表示某个条件对每个元素是否成立。

例：

1
a = [0 4 9 7 3 5] <= 4;

结果：

1
[1 1 0 0 1 0]

44. 创建逻辑向量

44.1 关系运算得到

1
a = [1 6 5] < 2;

结果：

1
[1 0 0]

44.2 直接指定

1
a = logical([1 0 0]);

如果不调用 logical，MATLAB 会创建 double 类型的 [1 0 0]。

45. Logical vector rules

例：

1
x = [5 9 2 4 3];
2
v = logical([1 0 1 0 1]);
3
xp = x(v);

结果：

1
xp = [5 2 3]

规则：

1. v 必须是 logical 类型。
2. 返回 x 中对应 v 为 true 的元素。
3. x 和 v 必须大小相同。

46. 逻辑向量例子

1
x = [1 4 3 6 2 8];
2
v = x > 3;

结果：

1
v = [0 1 0 1 0 1]

1
v = x < -3;

结果：

1
[0 0 0 0 0 0]

1
v = ~(~x);

由于 x 全非零，所以：

1
[1 1 1 1 1 1]

47. 逻辑向量应用：sin(x)

目标：把 sin(x) 中负值设为 0。

47.1 方法 1：for + if

1
x = 0:pi/1000:3*pi;
2
y = sin(x);
3

4
for jj = 1:length(y)
5
    if y(jj) < 0
6
        y(jj) = 0;
7
    end
8
end
9

10
plot(x,y,'k.')

47.2 方法 2：逻辑向量

1
x = 0:pi/1000:3*pi;
2
y = sin(x);
3

4
y = y.*(y > 0);
5

6
plot(x,y,'k.')

更高效。

47.3 绘图代码

1
x = 0:pi/100:4*pi;
2
y = sin(x);
3

4
yp = y.*(y>0);
5

6
pl = plot(x,y,'k-',x,yp,'r-');
7
set(gca,'FontSize',18)
8
set(pl(1),'LineWidth',10)
9
set(pl(2),'LineWidth',4)
10
axis([0 4*pi -1.2 1.2])
11
xlabel('x'), ylabel('y or y_p')
12
lg = legend('y=sin(x)',...
13
    'y_p=sin(x), y_p=0 for y<0');

48. 逻辑向量应用：tan(x)

原始代码：

1
x = 0*pi:pi/10000:pi;
2
y = tan(x);
3
plot(x,y)
4
xlim([0 pi])
5
xlabel('x')
6
ylabel('y')

由于 tan 在 pi/2 附近有渐近线，图像细节会被大值支配。

48.1 大值设为 0

1
x = 0*pi:pi/10000:pi;
2
y = tan(x);
3
y = y.*(abs(y)<15);
4
plot(x,y)

48.2 大值设为 NaN

1
x = 0*pi:pi/10000:pi;
2
y = tan(x);
3
y(abs(y)>15) = NaN;
4
plot(x,y)

NaN 不会被正常连线绘制，因此能断开渐近线附近的图。

48.3 删除大值

1
x = 0*pi:pi/10000:pi;
2
y = tan(x);
3
x = x(abs(y)<15);
4
y = y(abs(y)<15);
5
plot(x,y)

49. Example：Sieve of Eratosthenes 埃拉托色尼筛法

目标：找出 1 到 n 的所有质数。

质数定义：

• 只能被 1 和自身整除
• 必须是大于 1 的整数

算法：

1. 列出从 2 到 n 的数。
2. 当前列表第一个数一定是质数，把它加入质数列表。
3. 划掉该数的所有倍数。
4. 重复步骤 2–4，直到当前质数大于等于 sqrt(n)。
5. 剩余所有数都是质数。

该算法使用：

• logical expressions
• logical indexing

Lec 3：矩阵、多维数组、矩阵运算、稀疏矩阵、数组函数

1. What is a matrix?

矩阵是二维数组，由：

• rows 行
• columns 列

组成。

例如一个 3 行 3 列矩阵大小为：

1
3 x 3

2. Multidimensional Array / Matrix

数据类型可从低维到高维理解：

| 类型 | 含义 | | - | | | Scalar | 标量，一个数 | | Vector | 向量，一维数组 | | Matrix / Array | 矩阵 / 二维数组 | | Multidimensional Matrix / Array | 多维数组 |

三维数组可以理解为多个二维矩阵沿第三维堆叠。

3. Basic Array-related Operators

3.1 方括号

用于创建数组。

数组数据按行列组织，每行长度必须相同。

行内元素

用逗号或空格：

1
A = [1 2 3];
2
A = [1, 2, 3];

换行

用分号：

1
A = [1 2 3; 4 5 6];

结果是 2×3 矩阵。

3.2 冒号

用于创建大数组和等间距数据。

默认步长为 1

1
first_value:last_value

自定义步长

1
first_value:step:last_value

例：

1
[f:n]

生成：

1
[f, f+1, f+2, ..., k]

其中 k <= n 且 k+1 > n。 如果 f > n，则结果为空数组 []。

4. Creating arrays 创建数组

4.1 创建一维数组 / 向量

常间距向量

1
My_variable = [first:spacing:last];
2
A = [1:5];
3
A = [1:2:5];

也可以不写方括号：

1
My_variable = first:spacing:last;
2
x = 1:5;
3
x = 1:0.5:5;

linspace

1
My_variable = linspace(first, last, num);
2
B = linspace(1,50,10);

表示从 1 到 50 均匀取 10 个点。

任意向量

1
A = [1 3 0 -2 5 6 10];

5. Creating matrices 创建二维矩阵

5.1 手动创建

1
A = [1 2 3; 4 5 6];

结果：

1
1 2 3
2
4 5 6

5.2 内置函数

1
A = ones(2,5);
2
A = zeros(2,3);
3
A = eye(3);

6. Transpose operator 转置 '

1
A'

作用：

• 行变列
• 列变行
• m×n 矩阵变 n×m 矩阵
• 列向量变行向量
• 行向量变列向量

例：

1
A = [1 2 3; 4 5 6];
2
A'

结果：

1
1 4
2
2 5
3
3 6

7. 使用表达式初始化数组

数组初始化表达式可以包含：

• 算术运算
• 已定义数组的全部或部分

例：

1
a = [0 13*2];

结果：

1
a = [0 26]

1
b = [a(2) 13 a];

结果：

1
b = [26 13 0 26]

8. Accessing matrix indices 矩阵索引

8.1 二维下标访问

格式：

1
A(row, column)

例如：

1
A(3,1)

表示第 3 行第 1 列。

9. 单下标访问矩阵

矩阵也可以用单个下标访问。 MATLAB 按列展开矩阵，即 column-major order：

从左上角开始，先沿第一列向下，再第二列向下。

对 3×3 矩阵，单下标顺序：

1
1 4 7
2
2 5 8
3
3 6 9

也就是说 A(1) 是第一行第一列，A(2) 是第二行第一列。

10. Subarrays 子数组

10.1 向量

1
va(:)

所有元素。

1
va(m:n)

第 m 到第 n 个元素。

10.2 矩阵

10.3 end

end 表示对应维度的最大索引。

例：

1
a = [1 2 3 4 5 6 7 8 9];
2
a(7:end)

结果：

1
[7 8 9]

1
a(end)

结果：

1
9

10.4 同一表达式中 end 可代表不同维度

1
b = [1 2 3 4; 5 6 7 8; 9 10 11 12];
2
b(2:2:end,2:end)

含义：

• 行：从第 2 行到最后一行，每隔 2 行取一次
• 列：第 2 列到最后一列

结果：

1
[6 7 8]

11. Manipulating / creating arrays

包括：

• 添加元素 expanding
• 删除元素
• 内置数组函数

12. Resizing arrays 调整数组大小

12.1 length

返回元素个数，或数组最大维度长度。

1
length(v)

12.2 size(A)

返回矩阵大小：

1
[m,n] = size(A);

或者：

1
size(A)

返回：

1
[m n]

12.3 reshape(A,m,n)

重排 A，使其有 m 行 n 列。

1
B = reshape(A,m,n);

注意：元素总数必须一致。

12.4 diag(v)

如果输入是向量，创建以该向量为主对角线的方阵。

1
diag([1 2 3])

结果：

1
1 0 0
2
0 2 0
3
0 0 3

12.5 diag(A)

如果输入是矩阵，返回主对角线元素向量。

12.6 自动扩展

1
A(end+1) = 5;

自动添加一个元素。

12.7 删除行或列

1
A(:,2) = [];

删除第 2 列。

12.8 拼接数组

1
A = [B C; D E];

含义：

1
左上 B，右上 C
2
左下 D，右下 E

前提是维度能对齐。

13. 对角线元素

如何指定矩阵 A 的对角线元素？

方法 1：

1
a = [A(1,1) A(2,2) A(3,3)];

方法 2：单下标

1
a = A(1:5:9);

对 3×3 矩阵，主对角线单下标为 1,5,9。

方法 3：

1
a = diag(A);

工程中对角线元素常很重要，例如 stress tensor 中的 normal stresses。

14. Strings 字符串

MATLAB 字符串可以是单引号中的字符数组。

例：

1
'ad ef'
2
'CSCSI0040'
3
'Ruth Simmons'
4
'99.99% pure'
5
'Exclaiming text!'

14.1 字符串中包含单引号

用两个单引号表示一个单引号：

1
'text with a quote''in it'

14.2 char 函数

1
My_string = char('string 1','string 2','string 3');

创建字符数组。

14.3 字符串相关函数

例如：

1
lower
2
isspace
3
isletter

15. Matrix Math 矩阵数学

内容包括：

• addition / subtraction
• multiplication
• division
• array math
• array operations

16. Array addition and subtraction

规则：

1. 标量可以加到任意数组，标量会加到每个元素。
2. 数组之间加减必须维度相同，对应位置逐元素加减。

16.1 标量扩展 scalar expansion

1
A = [1 2; 3 4];
2
B = A + 2;

结果：

1
B = [3 4; 5 6]

16.2 数组维度必须一致

1
B = A + [3 4];

报错：

1
matrix dimensions must agree

1
B = A + [3 4; 5 6];

结果：

1
B = [4 6; 8 10]

17. Element-by-element operations

数组大小必须一致。

内置函数如：

1
sin(A)
2
log(A)

会对数组每个元素分别计算。

18. Matrix / array multiplication

18.1 标量乘数组

1
B = 2*A;

MATLAB 把标量扩展到数组每个元素。

18.2 矩阵乘法

如果：

1
A: 4 x 3
2
B: 3 x 2

则：

1
A*B: 4 x 2

要求前一个矩阵列数等于后一个矩阵行数。

18.3 * 与 .*

1
A = [1 2; 3 4];

矩阵乘法：

1
B = A*[1 2; 3 4];

结果：

1
B = [7 10; 15 22]

逐元素乘法：

1
B = A.*[1 2; 3 4];

结果：

1
B = [1 4; 9 16]

19. Vector multiplication 向量乘法

19.1 行向量 × 列向量

结果是标量，即 dot product。

19.2 列向量 × 行向量

结果是矩阵。

19.3 Cross product

叉乘是特殊操作，使用：

1
cross(a,b)

20. Matrix division 矩阵除法

20.1 Identity matrix 单位矩阵

单位矩阵 I 是 n×n 方阵：

• 对角线为 1
• 其他位置为 0

满足：

1
AI = IA = A

20.2 Inverse matrix 逆矩阵

如果：

1
BA = AB = I

则 B 是 A 的逆矩阵。

MATLAB：

1
inv(A)
2
A^-1

注意：

• A 必须是方阵
• A 必须可逆

21. Left division 左除

求：

1
AX = B

理论：

1
X = A^(-1)B

MATLAB 推荐：

1
X = A\B;

不要优先使用：

1
X = inv(A)*B;

因为 A\B 更快更稳定。

22. Right division 右除

求：

1
XC = D

理论：

1
X = D C^(-1)

MATLAB：

1
X = D/C;

文件：

1
L9.m

23. Gaussian Elimination 高斯消元

可用于解线性方程组：

1
AX = B

或：

1
XC = D

文件：

1
L9.m

24. Transposing arrays

转置是沿主对角线反射矩阵。

1
A'

行向量转置为列向量，列向量转置为行向量。

25. Rotating and flipping arrays

MATLAB 可以旋转和翻转矩阵。

25.1 图像和矩阵的关系

图像的垂直轴通常和用户定义矩阵的垂直轴方向相反。

因此处理图像和矩阵混合问题时，可能需要翻转。

25.2 生成图像示例代码

1
figure(1)
2
load('spine','X','map')
3
imagesc(X)
4
colormap(map)
5
set(gca,'FontSize',20)
6
set(gcf,'Color','w')
7
xlabel('Index')
8
ylabel('Index')
9

10
figure(2)
11
rndArr = rand(367,490);
12
surface(rndArr)
13
set(gca,'FontSize',20)
14
set(gcf,'Color','w')
15
xlabel('Index')
16
ylabel('Index')
17
axis tight
18
colormap('gray')
19
shading interp

26. flipdim

rot90、fliplr、flipud 要求二维矩阵。 flipdim 可以处理任意维度数组。

1
flipdim(A,1)
2
flipdim(A,2)
3
flipdim(A,3)

分别沿第 1、2、3 维翻转。

注：较新 MATLAB 中 flipdim 已逐渐被 flip(A,dim) 替代。

27. Built-in matrices 内置矩阵

MATLAB 有很多预定义矩阵函数。

27.1 指定尺寸

单个参数：

1
zeros(3)

创建 3×3 矩阵。

两个参数：

1
ones(2,3)

创建 2×3 矩阵。

某些函数支持 N 维数组：

1
rand(m,n,p)
2
randi(m,n,p)

27.2 常见函数

1
zeros(3)
2
ones(2,3)
3
magic(3)
4
eye(3)
5
rand(2,3,5)
6
randi(10,3,5,4)

28. Sparse matrices 稀疏矩阵

稀疏矩阵适用于大部分元素为 0 的矩阵。

应用：

• search algorithms
• atomistics
• stress / strain mapping
• FEM

优点：

• 减少内存
• 加快计算

28.1 常规创建大矩阵

假设矩阵大小 1000×2000，共 2 million 数据点。

非零元素：

1
(3,4) = 15
2
(100,1500) = 5
3
(1000,2000) = 9

常规方法：

1
A = zeros(1000,2000);
2
A(3,4) = 15;
3
A(100,1500) = 5;
4
A(1000,2000) = 9;

28.2 sparse 创建

1
A = sparse([3,100,1000],...
2
           [4,1500,2000],...
3
           [15,5,9],1000,2000);

语法：

1
sparse(row_indices, column_indices, values, m, n)

28.3 full 转 sparse

1
A = sparse(A);

28.4 节省效果

| Matrix type | Memory | Time to compute A.^2 | | — | -: | : | | Full | 15.2 MB | 4 ms | | Sparse | 0.015 MB | 0.04 ms |

29. Array size

常用函数：