Lecture 2
θ
和它的导数符号是通过 Julia 中的变量命名方式实现的
-
变量
θ
的输入:- 在 Julia 中,
θ
是一个合法的变量名,就像普通的字母x
或y
一样。 - 要输入
θ
,可以使用以下方法:- 在 Jupyter Notebook 或 Julia REPL 中直接键入
\theta
,然后按 Tab 键,Julia 会自动将其转化为θ
。 - 这是 Julia 支持的 Unicode 字符的一部分,可以直接作为变量名。
- 在 Jupyter Notebook 或 Julia REPL 中直接键入
- 在 Julia 中,
-
变量
θ̇
的输入:θ̇
是带点的变量(表示导数),它也是 Julia 的一个合法变量名。- 键入
\theta
,后按 Tab 键,键入\dot
,后按 Tab 键
-
变量
θ̈
的输入:θ̈
是加速度(带两个点的变量名)。- 键入
\theta
,后按 Tab 键,键入\ddot
,后按 Tab 键
θ
和θ̇
是 Julia 中的 Unicode 变量名,可以通过快捷输入实现。- 在代码中,这些符号只是普通变量名,它们的含义通过上下文和赋值来表达:
θ
是角度。θ̇
是角速度。θ̈
是角加速度,计算公式是 − ( g / l ) sin ( θ ) -(g/l) \sin(θ) −(g/l)sin(θ)。
.
使得操作符能够逐个元素地应用于整个数组或集合,而不需要显式地使用循环
在 Julia 中,等号前面的 .
是一个 广播(broadcasting) 操作符,它的作用是让一个操作应用于整个数组或容器,而不仅仅是单一的元素。这是 Julia 中非常常见和强大的特性。
x_hist[:,1] .= x0
x_hist[:,1]
:这是对x_hist
数组的访问,表示提取x_hist
的第一列(所有行,第一列)。.=
:这是广播赋值操作符。它表示对x_hist[:,1]
的每个元素执行赋值操作,等价于逐元素地把x0
的每个值赋给x_hist[:,1]
的每个元素。x0
:这是一个向量,通常是初始状态向量,比如[0.1; 0]
,表示摆的初始角度和角速度。
示例
假设 x0
是 [0.1, 0]
,而 x_hist[:,1]
是一个长度为 2 的列向量(例如 [0, 0]
),那么:
-
无广播(普通赋值):
如果直接使用=
(没有点操作符),会发生维度不匹配的错误,因为你不能将一个列向量直接赋值给数组的某一列。x_hist[:,1] = x0 # 这会报错
-
使用广播赋值:
使用.=
后,x0
的每个元素都会依次赋值给x_hist[:,1]
的每个元素:x_hist[:,1] .= x0 # 正确:x0 的每个元素逐个赋值给 x_hist[:,1]
总结
.
操作符使得赋值操作可以应用到整个数组或容器的每个元素,而不仅仅是简单的元素赋值。.=
用于广播赋值,使得多个值可以逐元素地赋给目标数组或集合。
这种广播机制使得在 Julia 中处理数组和矩阵时更加简洁高效。
解释匿名函数 x -> fd_pendulum_rk4(x, 0.1)
在 Julia 中,匿名函数是一种无需命名的函数,通常用于临时计算。形式为:
x -> expression
其中,x
是输入参数,expression
是处理该参数的表达式。
代码部分
x -> fd_pendulum_rk4(x, 0.1)
这是一个匿名函数,具体含义如下:
-
输入参数:
x
是匿名函数的输入参数,代表当前状态变量,通常为一个向量。例如,对于摆动系统,x
可能是 x = [ angle ; angular velocity ] x = [\text{angle}; \text{angular velocity}] x=[angle;angular velocity](角度和角速度)。
-
调用
fd_pendulum_rk4
:- 匿名函数内部调用了
fd_pendulum_rk4
函数,这是一个实现四阶龙格-库塔(RK4)积分方法的函数。 fd_pendulum_rk4(x, 0.1)
的两个参数:x
: 当前状态(例如 x = [ 角度 ; 角速度 ] x = [\text{角度}; \text{角速度}] x=[角度;角速度])。0.1
: 时间步长 h h h,表示模拟的离散时间间隔。
- 匿名函数内部调用了
-
输出结果:
- 函数返回 RK4 方法计算的状态更新结果,即从状态 x x x 经一步积分后的新状态 x n + 1 x_{n+1} xn+1。
- 该结果是一个新的状态向量,表示在给定时间步长下,系统从状态 x x x 演化到的下一个状态。
用途
匿名函数 x -> fd_pendulum_rk4(x, 0.1)
的核心作用是:
- 将输入状态 x x x(如摆的当前角度和角速度)映射为通过 RK4 方法计算得到的下一个状态。
- 在调用
ForwardDiff.jacobian
时,匿名函数为ForwardDiff
提供了所需的输入-输出关系。
示例
假设:
- 初始状态为 x = [ π / 4 ; 0 ] x = [\pi/4; 0] x=[π/4;0](角度 π / 4 \pi/4 π/4,角速度 0)。
- 时间步长 h = 0.1 h = 0.1 h=0.1。
匿名函数的计算流程如下:
# 定义匿名函数
f = x -> fd_pendulum_rk4(x, 0.1)# 调用匿名函数
new_state = f([pi/4; 0])
结果 new_state
是通过 RK4 积分计算出的下一步状态(新的角度和角速度)。
为什么使用匿名函数?
-
简洁性:
- 不需要显式定义一个新函数,而是直接将
fd_pendulum_rk4
封装成满足特定需求的函数(固定步长为 0.1)。
- 不需要显式定义一个新函数,而是直接将
-
灵活性:
- 匿名函数可以动态封装不同的参数和逻辑。例如,步长可以通过匿名函数灵活指定。
-
与
ForwardDiff.jacobian
配合:ForwardDiff.jacobian
需要输入一个函数,该函数的输入是状态 x x x,输出是对应的更新结果。匿名函数很好地满足这一要求。
总结
- 匿名函数
x -> fd_pendulum_rk4(x, 0.1)
将状态变量 x x x 映射为通过 RK4 方法计算得到的下一个状态。 - 它的主要作用是为
ForwardDiff.jacobian
提供输入-输出映射关系,以计算状态更新过程的雅可比矩阵。
Lecture 3
∂r \ r
表示 矩阵左除
在 Julia 中,表达式 ∂r \ r
表示 矩阵左除,也就是 求解线性方程组 的一种简洁方式。
线性方程组的求解
线性方程组的一般形式是:
A ⋅ x = b A \cdot x = b A⋅x=b
其中:
- A A A 是系数矩阵(这里对应
∂r
)。 - x x x 是未知量向量(这里对应 Δ x \Delta x Δx)。
- b b b 是右侧的已知向量(这里对应
r
)。
A \ x
的含义是 求解 x x x 的值,即:
x = A − 1 ⋅ b x = A^{-1} \cdot b x=A−1⋅b
- 这相当于将 A A A 的逆矩阵 A − 1 A^{-1} A−1 左乘到 b b b 上,求解 x x x 的值。
但是,显式求逆(即 A − 1 A^{-1} A−1)的计算代价很高,且可能会引入数值不稳定性。因此,∂r \ r
使用了一种数值高效的方式解决这个问题。
Julia 的 \
运算符
- 在 Julia 中,
A \ b
是求解线性方程组 A ⋅ x = b A \cdot x = b A⋅x=b 的符号,表示“将矩阵 A A A 左除向量 b b b”。 - 实现时,Julia 使用优化的数值线性代数方法(如 LU 分解、QR 分解或 Cholesky 分解)来高效求解,而不是直接计算矩阵的逆。
数值计算的优势
- 高效性: 求解 A ⋅ x = b A \cdot x = b A⋅x=b 的方法通常比显式逆矩阵的计算更高效。
- 数值稳定性: 显式计算逆矩阵可能会导致数值不稳定性(尤其当矩阵接近奇异时),而直接求解方程组能够减少误差。
- 灵活性: Julia 的
\
运算符会自动选择最适合的分解算法(如 LU、QR 或其他方法)来解决问题,适用于稠密矩阵或稀疏矩阵。
梯度和雅可比矩阵在 Julia 中的使用规则
- 梯度 (
gradient
): 用于标量值函数,返回一个列向量。 - 雅可比矩阵 (
jacobian
): 用于矢量值函数,返回一个矩阵。 - 在实际使用中,必须明确函数的输入和输出维度,误用可能导致报错。
1. 梯度(gradient)
梯度是用于标量值函数的,它返回的是一个列向量。
例子
using ForwardDiff# 定义一个标量函数
f(x) = x[1]^2 + x[2]^2 + x[3]^2# 对 f 求梯度
x = [1.0, 2.0, 3.0]
grad = ForwardDiff.gradient(f, x)
println("梯度: ", grad)
输出
- 解析:
- 函数 f ( x ) = x 1 2 + x 2 2 + x 3 2 f(x) = x_1^2 + x_2^2 + x_3^2 f(x)=x12+x22+x32。
- 梯度为 ∇ f = [ ∂ f ∂ x 1 , ∂ f ∂ x 2 , ∂ f ∂ x 3 ] T = [ 2 x 1 , 2 x 2 , 2 x 3 ] T \nabla f = [\frac{\partial f}{\partial x_1}, \frac{\partial f}{\partial x_2}, \frac{\partial f}{\partial x_3}]^T = [2x_1, 2x_2, 2x_3]^T ∇f=[∂x1∂f,∂x2∂f,∂x3∂f]T=[2x1,2x2,2x3]T。
- 结果为
[2.0, 4.0, 6.0]
。
2. 雅可比矩阵(Jacobian)
雅可比矩阵用于矢量值函数,返回的是一个矩阵。
例子
using ForwardDiff# 定义一个矢量值函数
g(x) = [x[1]^2, x[1]*x[2], x[2]^2]# 对 g 求雅可比矩阵
x = [1.0, 2.0]
jacobian = ForwardDiff.jacobian(g, x)
println("雅可比矩阵: ")
println(jacobian)
输出
- 解析:
-
函数 g ( x ) = [ x 1 2 , x 1 x 2 , x 2 2 ] g(x) = [x_1^2, x_1 x_2, x_2^2] g(x)=[x12,x1x2,x22]。
-
雅可比矩阵为:
J = [ ∂ g 1 ∂ x 1 ∂ g 1 ∂ x 2 ∂ g 2 ∂ x 1 ∂ g 2 ∂ x 2 ∂ g 3 ∂ x 1 ∂ g 3 ∂ x 2 ] J = \begin{bmatrix} \frac{\partial g_1}{\partial x_1} & \frac{\partial g_1}{\partial x_2} \\ \frac{\partial g_2}{\partial x_1} & \frac{\partial g_2}{\partial x_2} \\ \frac{\partial g_3}{\partial x_1} & \frac{\partial g_3}{\partial x_2} \end{bmatrix} J= ∂x1∂g1∂x1∂g2∂x1∂g3∂x2∂g1∂x2∂g2∂x2∂g3 = [ 2 x 1 0 x 2 x 1 0 2 x 2 ] \begin{bmatrix} 2x_1 & 0 \\ x_2 & x_1 \\ 0 & 2x_2 \end{bmatrix} 2x1x200x12x2
-
在 x = [ 1.0 , 2.0 ] x = [1.0, 2.0] x=[1.0,2.0] 时,结果为:
[ 2.0 0.0 2.0 1.0 0.0 4.0 ] \begin{bmatrix} 2.0 & 0.0 \\ 2.0 & 1.0 \\ 0.0 & 4.0 \end{bmatrix} 2.02.00.00.01.04.0
-
3. 错误调用的情况
错误调用 gradient
对矢量值函数
如果尝试对矢量值函数调用 gradient
会导致报错,因为梯度只适用于标量值函数。
例子
g(x) = [x[1]^2, x[1]*x[2], x[2]^2]
x = [1.0, 2.0]
grad = ForwardDiff.gradient(g, x) # 错误
报错信息
- 原因:
gradient
只能对标量值函数使用,而这里的 g ( x ) g(x) g(x) 是矢量值函数。
错误调用 jacobian
对标量值函数
如果尝试对标量值函数调用 jacobian
,理论上应返回梯度的转置,但通常会导致报错。
例子
f(x) = x[1]^2 + x[2]^2
x = [1.0, 2.0]
jacobian = ForwardDiff.jacobian(f, x) # 错误
报错信息
- 原因:
jacobian
期望输入是矢量值函数,而这里的 f ( x ) f(x) f(x) 是标量值函数。
Lecture 4
Kronecker 积(Kronecker Product)
Kronecker 积的定义
给定两个矩阵 A ∈ R m × n A \in \mathbb{R}^{m \times n} A∈Rm×n 和 B ∈ R p × q B \in \mathbb{R}^{p \times q} B∈Rp×q,它们的 Kronecker 积 A ⊗ B A \otimes B A⊗B 是一个大小为 ( m p ) × ( n q ) (mp) \times (nq) (mp)×(nq) 的矩阵。具体构造规则如下:
- A ⊗ B A \otimes B A⊗B 将矩阵 A A A 的每个元素 a i j a_{ij} aij 替换为该元素与矩阵 B B B 的乘积 a i j B a_{ij}B aijB。
数学表达为:
A ⊗ B = [ a 11 B a 12 B … a 1 n B a 21 B a 22 B … a 2 n B ⋮ ⋮ ⋮ a m 1 B a m 2 B … a m n B ] A \otimes B = \begin{bmatrix} a_{11} B & a_{12} B & \dots & a_{1n} B \\ a_{21} B & a_{22} B & \dots & a_{2n} B \\ \vdots & \vdots & & \vdots \\ a_{m1} B & a_{m2} B & \dots & a_{mn} B \\ \end{bmatrix} A⊗B= a11Ba21B⋮am1Ba12Ba22B⋮am2B………a1nBa2nB⋮amnB
Kronecker 积的计算方法
假设:
A = [ 1 2 3 4 ] , B = [ 0 5 6 7 ] A = \begin{bmatrix} 1 & 2 \\ 3 & 4 \\ \end{bmatrix}, \quad B = \begin{bmatrix} 0 & 5 \\ 6 & 7 \\ \end{bmatrix} A=[1324],B=[0657]
计算 A ⊗ B A \otimes B A⊗B:
A ⊗ B = [ 1 ⋅ B 2 ⋅ B 3 ⋅ B 4 ⋅ B ] = [ [ 0 5 6 7 ] [ 0 10 12 14 ] [ 0 15 18 21 ] [ 0 20 24 28 ] ] A \otimes B = \begin{bmatrix} 1 \cdot B & 2 \cdot B \\ 3 \cdot B & 4 \cdot B \\ \end{bmatrix}= \begin{bmatrix} \begin{bmatrix} 0 & 5 \\ 6 & 7 \\ \end{bmatrix} & \begin{bmatrix} 0 & 10 \\ 12 & 14 \\ \end{bmatrix} \\ \begin{bmatrix} 0 & 15 \\ 18 & 21 \\ \end{bmatrix} & \begin{bmatrix} 0 & 20 \\ 24 & 28 \\ \end{bmatrix} \\ \end{bmatrix} A⊗B=[1⋅B3⋅B2⋅B4⋅B]= [0657][0181521][0121014][0242028]
最终结果为:
A ⊗ B = [ 0 5 0 10 6 7 12 14 0 15 0 20 18 21 24 28 ] A \otimes B = \begin{bmatrix} 0 & 5 & 0 & 10 \\ 6 & 7 & 12 & 14 \\ 0 & 15 & 0 & 20 \\ 18 & 21 & 24 & 28 \\ \end{bmatrix} A⊗B= 0601857152101202410142028
Kronecker 积的性质
- 尺寸:
如果 A A A 是 m × n m \times n m×n, B B B 是 p × q p \times q p×q,那么 A ⊗ B A \otimes B A⊗B 的大小为 ( m p ) × ( n q ) (mp) \times (nq) (mp)×(nq)。 - 分布律:
( A + C ) ⊗ B = A ⊗ B + C ⊗ B (A + C) \otimes B = A \otimes B + C \otimes B (A+C)⊗B=A⊗B+C⊗B - 结合律:
( A ⊗ B ) ⊗ C = A ⊗ ( B ⊗ C ) (A \otimes B) \otimes C = A \otimes (B \otimes C) (A⊗B)⊗C=A⊗(B⊗C) - 与标量的关系:
( α A ) ⊗ B = A ⊗ ( α B ) (\alpha A) \otimes B = A \otimes (\alpha B) (αA)⊗B=A⊗(αB)
Kronecker 积的应用
-
生成重复矩阵:
- 在 Julia 中,
kron(ones(m), A)
会生成一个矩阵,其中矩阵 A A A 的每一行重复 m m m 次。 - 类似地,
kron(ones(m)', A)
会生成一个矩阵,其中矩阵 A A A 的每一列重复 m m m 次。
- 在 Julia 中,
-
向量化操作:
- Kronecker 积常用于将向量化表达与矩阵展开结合。比如,将矩阵 A A A 的每一项与另一个矩阵 B B B 关联。
-
量子计算:
- Kronecker 积在量子力学中用于描述复合量子系统的状态和操作,比如计算张量积态。
-
系统理论和信号处理:
- Kronecker 积用于构造大规模系统矩阵,特别是在多维信号处理中的应用。
在 Julia 中的实现
在 Julia 中,kron
函数用于计算 Kronecker 积,语法为:
C = kron(A, B)
在equality-constraints.ipynb代码中的作用
在 plot_landscape
函数中,Kronecker 积被用来生成网格点:
kron(ones(Nsamp), LinRange(-4, 4, Nsamp)')
:- 生成一个矩阵,其中每一行是从 − 4 -4 −4 到 4 4 4 的序列,表示 x x x 坐标。
kron(ones(Nsamp)', LinRange(-4, 4, Nsamp))
:- 生成一个矩阵,其中每一列是从 − 4 -4 −4 到 4 4 4 的序列,表示 y y y 坐标。
这样,通过 Kronecker 积可以快速构造二维网格,用于绘制等高线图。