目前市面上的Julia教程很多,但大部分教程都是从0开始,对我这种有其他语言编程基础的十分不友好,我真的不关注你的Unicode支持有多棒好吗,类似这种东西在教程里大讲特讲对我来说属实没必要。另一部分则是过于高深,各种高级特性满天飞,本想着先拿来写点东西,看到奇怪的语法当场被劝退。恰好最近在学数值分析,借这个机会刚好练习一下Julia的使用。我的使用仅关注功能的实现,在语法上仅限能跑通,对其他细节则不去关注。
本文的代码可以在这里看到
AppliedNumericalAnalysisgithub.com
先介绍一下本文用到的Julia语法
- Julia的索引从1开始,与Python、C++不一样
- Array即我们需要的矩阵
x = [1. 2. 3.; 2 7 5; 1 4 9] # 3 * 3的矩阵
b = [1; 6; -3.] # 长度为3的列向量
3. size用来获取矩阵的大小
size(x) # 返回二维矩阵的行数和列数
size(x, 1) # 返回矩阵的行数
size(x, 2) # 返回矩阵的列数
4. 矩阵的索引方法和matlab、python类似
x = [1 2 3; 4 5 6; 7 8 9]
x[1, 1] # 1
x[1, :] # 取第一行的值,返回一个列向量
x[end, end] # 9
x[:, 2] # 取第二列的值,返回一个列向量
5. 用.*可以实现对应元素相乘(和matlab一样)
6. Julia的map和python类似,第一个参数为表达式,第二个参数为对象
高斯消去法
高斯消去法是用来求解线性方程组的经典解法,其基本思路是对增广矩阵做行变换化为上三角矩阵后再回代求得方程组的解。以下均假设方程组有解。
首先定义函数的签名
function gauss_elimination(X::Array, b::Array)
其中,X是n阶方阵,b是n维列向量。
为了便于运算,我们将方阵和列向量合并为增广矩阵
A = hcat(x, b)
row_num = size(A, 1)
接下来进行消去运算
for i in 1:row_num-1
pivot = A[i, i] # 对角线上元素
for j in i+1:row_num
base = A[j, i] / pivot # 确定消去该行第一个非零元素的系数
A[j,:] = A[j,:] - (base .* A[i,:]) # 对该行所有元素减去系数乘第i行的对应值
end
end
以如下线性方程组为例
上述代码模拟的即是消去的过程
在得到变换后的增广矩阵后即可进行回代
b = A[:, end]
b[end] /= A[end, end-1] # 最后一行直接计算即可
for i in row_num-1:-1:1 # 对其他行,先将除第一个非0列以外的其他列归零
pivot = A[i,i]
b[i] -= sum(A[i,i+1:end-1] .* b[i+1:end])
b[i] /= pivot
end
为了遵循课本的省内存的需求,将b作为结果返回。
高斯主元素法
高斯主元素法,又名列主元消去法,主要是为了解决在工程计算中由于计算机精度不足时除数过小所造成的精度丢失问题。其主要思想是,在每次消元时选取待选元素中最大的值作为当前目标。具体的思路可以在这里看到。
主元素法的基本思路和普通的消去法一致,仅多了选择最大列主元一项
target = map(item -> abs(item), A[i:end, i]) # 得到待选列
max_val = maximum(target) # 待选列最大值
max_index = findfirst(target .== max_val) # 待选列索引
if max_index != 1 # 若当前行的列值不是最大,则交换两行
A[i, :], A[max_index, :] = A[max_index, :], A[i, :]
end