Julia使用nlsolve求解非线性方程组时遇MethodError问题求助

问题重现 #

运行代码调用sys_tHE_AB!函数时触发如下错误：

ERROR: MethodError: no method matching +(::Float64, ::Array{Float64, 3})
For element-wise addition, use broadcasting with dot syntax: scalar .+ array
Closest candidates are:
  +(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
  +(::T, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383
  +(::Union{Float16, Float32, Float64}, ::BigFloat) at mpfr.jl:414

错误根源 #

1. 维度不匹配：主代码中直接将三维数组k_e、二维矩阵fe/Fe传入sys_tHE_AB!，但nlsolve求解的是4个标量未知数的方程组，要求F的每个元素必须是标量等式。标量与多维数组直接进行加减操作，触发类型不匹配错误。
2. 多余代码：函数末尾的H = [F[1] F[2] F[3] F[4]]无意义，nlsolve仅需要修改传入的F数组，不需要返回额外变量。

修复步骤 #

1. 主代码：传入对应索引的标量参数 #

在循环中，根据当前的c、d、j索引，提取k_e、fe、Fe的标量值传入函数，而非整个数组/矩阵：

# 求解阈值与价值函数参数
for j in 1:3
   for c in 1:N
      for d in 1:N+1
         if d != c
             # 省略部分初始化代码
                # 传入对应循环索引的标量，而非整个多维数组
                f_thresholds_HE!(F, Z_AB) = sys_tHE_AB!(F, Z_AB, k_e[c,d,j], fe[c,d], Fe[c,d], param, j)

                initial_x = init_HE[i,j,k,:]
                res = nlsolve(f_thresholds_HE!, initial_x, method = :trust_region, xtol=10e-5, ftol=10e-8, iterations=10000, store_trace=true, extended_trace=true)
            end
        end
end

注：索引顺序k_e[c,d,j]需根据你的数组实际维度定义调整，确保取到的是单个标量值。

2. 修正sys_tHE_AB!函数 #

更新参数类型为标量，移除多余代码，确保所有运算为标量间操作：

# 此函数计算4方程系统，求解出口阈值及参数B、A
function sys_tHE_AB!(F::AbstractArray, Z_AB::AbstractArray, k_e::Float64, fe::Float64, Fe::Float64, param::Parameter, j::Int64)
    μ = param.μ
    σ = param.σ
    μj = μ[j]
    σj = σ[j]
    ρ = param.ρ
    δ = param.δ

    β =  sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    α = -sqrt(2*(ρ + δ)/(σj^2) + (μj/(σj^2) - 1/2)^2) - (μj/(σj^2) - 1/2)
    disc = ρ + δ - μj

    Y_OE = Z_AB[1]
    Y_EO = Z_AB[2]
    Ae   = Z_AB[3]
    Bo   = Z_AB[4]

    # 所有运算均为标量操作，无维度不匹配问题
    F[1] = Ae*Y_OE^α + k_e*Y_OE/disc - fe/disc - Fe - Bo*Y_OE^β
    F[2] = Ae*Y_EO^α + k_e*Y_EO/disc - fe/disc - Bo*Y_EO^β
    F[3] = α*Ae*Y_OE^(α-1) + k_e/disc - β*Bo*Y_OE^(β-1)
    F[4] = α*Ae*Y_EO^(α-1) + k_e/disc - β*Bo*Y_EO^(β-1)
end

额外说明 #

如果需要批量求解多组方程组，可保留多维数组传入，但必须对所有标量与数组的运算使用广播语法（运算符前加.），例如：

F[1] = Ae*Y_OE^α .+ k_e*Y_OE/disc .- fe/disc .- Fe .- Bo*Y_OE^β

同时需确保F的维度与运算结果匹配，适配nlsolve的批量求解逻辑。

内容的提问来源于stack exchange，提问作者Coco de Stael