专栏
标签
Jacobi迭代
一般问题
发布于 2024-12-06 14:50:31
查看 29过去569天

function GaussMatrix(A, B)
    aug_matrix = hcat(A, B)
    (m, n) = size(A)
    for i = 1:n - 1
        max_row = argmax(abs.(aug_matrix[i:m, i]), dims = 1)[1] + i - 1
        if max_row!= i
            aug_matrix[[i, max_row], :] = aug_matrix[[max_row, i], :]
        end
        pivot = aug_matrix[i, i]
        if pivot == 0
            println("不能用高斯消去")
            return zeros(Float64, n)
        end
        multiplier_matrix = Matrix{Float64}(I, m, m)
        for j = i + 1:m
            multiplier_matrix[j, i] = -aug_matrix[j, i] / pivot
        end
        aug_matrix = multiplier_matrix * aug_matrix
    end
    x = zeros(Float64, n)
    x[n] = aug_matrix[n, end] / aug_matrix[n, n]
    for k = (n - 1):-1:1
        s = dot(aug_matrix[k, (k + 1):n], x[(k + 1):n])
        x[k] = (aug_matrix[k, end] - s) / aug_matrix[k, k]
    end
    return x
end
function JacobiMatrix(A, B, tol = 1e-8)
    (m, n) = size(A)
    if m!= n
        error("系数矩阵A必须是方阵")
    end
    if size(B, 2)!= 1
        B = reshape(B, n, 1)
    elseif size(B, 1)!= n
        error("常数项向量B的维度与系数矩阵A的行数不匹配")
    end
    row_sums = sum(abs.(A), dims = 2) - diagm(diag(abs.(A)))
    x = zeros(Float64, n, 1)
    x = B./ (diag(A) + row_sums[:, 1])
    max_iterations = 1000
    iteration_count = 0
    while true
        x_new = zeros(Float64, n, 1)
        Ax = A * x[:, 1]
        diagAx = diagm(diag(A)) * x[:, 1]
        diff = B - (Ax - diagAx)
        x_new = (diff / diag(A))[:, newaxis]
        error = norm(x_new - x)
        if error < tol || iteration_count >= max_iterations
            return x_new[:, 1]
        end
        x = x_new
        iteration_count += 1
    end
end
A = [10 -1 2 0; -1 11 -1 3; 2 -1 10 -1; 0 3 -1 8]
B = [6; 25; -11; 15]
solution_gauss_matrix = GaussMatrix(A, B)
println("The solution of Gauss (matrix version) is")
for i in 1:length(solution_gauss_matrix)
    println("x$i = ", solution_gauss_matrix[i])
end
solution_jacobi_matrix = JacobiMatrix(A, B)
println("The solution of Jacobi (matrix version) is")
for i in 1:length(solution_jacobi_matrix)
    println("x$i = ", solution_jacobi_matrix[i])
end
Jacobi迭代Julia脚本报错ERROR: DimensionMismatch: dimensions must match: a has dims (Base.OneTo(4), Base.OneTo(1)), b has dims (Base.OneTo(4), Base.OneTo(4)), mismatch at 2,这个怎么修改代码?

所属专栏:Julia语言
产品信息:Syslab科学计算环境
科学计算

全部回答 1

发布于 2024-12-07 09:50:38

您好,您的第38行代码,减号左边的维度是(4,1),减号右边的维度是(4,4),维度不同,无法直接相减,需要使用.-的形式

用户
和原帖交流更多问题细节吧,去
我要发帖 我要发帖
资料中心 资料中心
查看更多>
热门帖子 热门帖子
主要贡献者 主要贡献者
过去7天