可以用三种方法构造该三对角块方阵吗?
1.稀疏对角阵库函数spdiagm
2.矩阵水平连接竖直连接
3.张量积函数kron
源代码如下
using TyBase, TyMath, TyPlot
J = [-4 1 0; 1 -4 1; 0 1 -4]
E = [1 0 0; 0 1 0; 0 0 1]
Y = [0 0 0; 0 0 0; 0 0 0]
println("方法1: 使用spdiagm构造")
n = 3
block_size = 3
total_size = n * block_size
N1 = zeros(total_size, total_size)
for i in 1:n
N1[(i-1)block_size+1:iblock_size, (i-1)block_size+1:iblock_size] = J
if i < n
N1[(i-1)block_size+1:iblock_size, i*block_size+1:(i+1)block_size] = E
N1[iblock_size+1:(i+1)*block_size, (i-1)block_size+1:iblock_size] = E
end
end
diag_blocks = [J, J, J]
upper_blocks = [E, E]
lower_blocks = [E, E]
N1_sparse = spdiagm(0 => vcat([diag_blocks[i][j,j] for i in 1:n for j in 1:block_size]...),
block_size => vcat([upper_blocks[i][j,j] for i in 1:n-1 for j in 1:block_size]...),
-block_size => vcat([lower_blocks[i][j,j] for i in 1:n-1 for j in 1:block_size]...))
display(N1)
println("\n方法2: 使用矩阵水平和竖直连接")
row1 = hcat(J, E, Y)
row2 = hcat(E, J, E)
row3 = hcat(Y, E, J)
N2 = vcat(row1, row2, row3)
display(N2)
println("\n方法3: 使用张量积kron")
T = [-1 1 0; 1 -1 1; 0 1 -1]
I3 = [1 0 0; 0 1 0; 0 0 1]
N3 = kron(I3, J) + kron(T, Y)
display(N3)
全部回答 3
我只会构造三对角矩阵:
Mworks代码:
clc(),clear()
n=5;vo=ones(n)2
v_1=-ones(n-1);v1=3ones(n-1)
A=spdiagm(0=>vo,-1=>v_1,1=>v1)
print("A=")
display(A)
A=full(A)
print("A=");display(A)
运行结果:
A=5×5 SparseMatrixCSC{Float64, Int64} with 13 stored entries:
2.0 3.0 ⋅ ⋅ ⋅
-1.0 2.0 3.0 ⋅ ⋅
⋅ -1.0 2.0 3.0 ⋅
⋅ ⋅ -1.0 2.0 3.0
⋅ ⋅ ⋅ -1.0 2.0
A=5×5 Matrix{Float64}:
2.0 3.0 0.0 0.0 0.0
-1.0 2.0 3.0 0.0 0.0
0.0 -1.0 2.0 3.0 0.0
0.0 0.0 -1.0 2.0 3.0
0.0 0.0 0.0 -1.0 2.0
源代码如下
using TyBase, TyMath, TyPlot
J = [-4 1 0; 1 -4 1; 0 1 -4]
E = [1 0 0; 0 1 0; 0 0 1]
Y = [0 0 0; 0 0 0; 0 0 0]
println("方法1: 使用spdiagm构造")
n = 3
block_size = 3
total_size = n * block_size
N1 = zeros(total_size, total_size)
for i in 1:n
N1[(i-1)block_size+1:iblock_size, (i-1)block_size+1:iblock_size] = J
if i < n
N1[(i-1)block_size+1:iblock_size, i*block_size+1:(i+1)block_size] = E
N1[iblock_size+1:(i+1)*block_size, (i-1)block_size+1:iblock_size] = E
end
end
diag_blocks = [J, J, J]
upper_blocks = [E, E]
lower_blocks = [E, E]
N1_sparse = spdiagm(0 => vcat([diag_blocks[i][j,j] for i in 1:n for j in 1:block_size]...),
block_size => vcat([upper_blocks[i][j,j] for i in 1:n-1 for j in 1:block_size]...),
-block_size => vcat([lower_blocks[i][j,j] for i in 1:n-1 for j in 1:block_size]...))
display(N1)
println("\n方法2: 使用矩阵水平和竖直连接")
row1 = hcat(J, E, Y)
row2 = hcat(E, J, E)
row3 = hcat(Y, E, J)
N2 = vcat(row1, row2, row3)
display(N2)
println("\n方法3: 使用张量积kron")
T = [-1 1 0; 1 -1 1; 0 1 -1]
I3 = [1 0 0; 0 1 0; 0 0 1]
N3 = kron(I3, J) + kron(T, Y)
display(N3)
运行结果:
julia> 正在运行 2.6.jl
方法1: 使用spdiagm构造
9×9 Matrix{Float64}:
-4.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0
1.0 -4.0 1.0 0.0 1.0 0.0 0.0 0.0 0.0
0.0 1.0 -4.0 0.0 0.0 1.0 0.0 0.0 0.0
1.0 0.0 0.0 -4.0 1.0 0.0 1.0 0.0 0.0
0.0 1.0 0.0 1.0 -4.0 1.0 0.0 1.0 0.0
0.0 0.0 1.0 0.0 1.0 -4.0 0.0 0.0 1.0
0.0 0.0 0.0 1.0 0.0 0.0 -4.0 1.0 0.0
0.0 0.0 0.0 0.0 1.0 0.0 1.0 -4.0 1.0
0.0 0.0 0.0 0.0 0.0 1.0 0.0 1.0 -4.0
方法2: 使用矩阵水平和竖直连接
9×9 Matrix{Int64}:
-4 1 0 1 0 0 0 0 0
1 -4 1 0 1 0 0 0 0
0 1 -4 0 0 1 0 0 0
1 0 0 -4 1 0 1 0 0
0 1 0 1 -4 1 0 1 0
0 0 1 0 1 -4 0 0 1
0 0 0 1 0 0 -4 1 0
0 0 0 0 1 0 1 -4 1
0 0 0 0 0 1 0 1 -4
方法3: 使用张量积kron
9×9 Matrix{Int64}:
-4 1 0 0 0 0 0 0 0
1 -4 1 0 0 0 0 0 0
0 1 -4 0 0 0 0 0 0
0 0 0 -4 1 0 0 0 0
0 0 0 1 -4 1 0 0 0
0 0 0 0 1 -4 0 0 0
0 0 0 0 0 0 -4 1 0
0 0 0 0 0 0 1 -4 1
0 0 0 0 0 0 0 1 -4