首页 > 解决方案 > 非线性求解器总是产生零残差

问题描述

我正在学习如何解决非线性方程和 Julia (1.3.1),并想问我们应该如何使用 NLsolve。

作为第一步,我尝试以下操作:

using NLsolve
uni = 1
z = 3
x = [3,2,4,5,6]
y = z .+ x
print(y)
function g!(F, x)
    F[1:uni+4] = y .- z - x
end
nlsolve(g!, [0.5,0.9,1,2,3])

而且,我确认它的工作原理如下:

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 0.9, 1.0, 2.0, 3.0]
 * Zero: [2.999999999996542, 2.000000000003876, 4.000000000008193, 4.999999999990685, 5.999999999990221]
 * Inf-norm of residuals: 0.000000
 * Iterations: 2
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 3
 * Jacobian Calls (df/dx): 3

然后,我尝试了一个更复杂的模型,如下所示

using SpecialFunctions, NLsolve, Random
Random.seed!(1234)
S = 2

# Setting parameters
ε = 1.3 
β = 0.4 
γ = gamma((ε-1)/ε) 
T = rand(S) 
E = rand(S)
B = rand(S)
w = rand(S)
Q = rand(S)
d = [1 2 ;2 1 ]

# Construct a model

rvector = T.*Q.^(1-β).*B.^ε
svector = E.* w.^ε


Φ_all = (sum(sum(rvector * svector' .* d )))
π = rvector * svector' .* d ./ Φ_all 

# These two are outcome the model 
πR = (sum(π,dims=1))' 
πM = sum(π,dims=2) 

# E is now set as unknown and we want to estimate it given the outcome of the model

function f!(Res, Unknown)


    rvector = T.*Q.^(1-β).*B.^ε
    svector = Unknown[1:S].* w.^ε 


    Φ_all = (sum(sum(rvector * svector' .* d )))
    π = rvector * svector' .* d ./ Φ_all
    Res = ones(S+S)
    Res[1:S] = πR - (sum(π,dims=1))' 
    Res[S+1:S+S] = πM - sum(π,dims=2) 

end


nlsolve(f!, [0.5,0.6])

这段代码产生了一个奇怪的结果,如下所示。

Results of Nonlinear Solver Algorithm
 * Algorithm: Trust-region with dogleg and autoscaling
 * Starting Point: [0.5, 0.6]
 * Zero: [0.5, 0.6]
 * Inf-norm of residuals: 0.000000
 * Iterations: 0
 * Convergence: true
   * |x - x'| < 0.0e+00: false
   * |f(x)| < 1.0e-08: true
 * Function Calls (f): 1
 * Jacobian Calls (df/dx): 1

因此,本质上,该函数总是产生 0 作为返回值,因此初始输入总是成为解决方案。我无法理解为什么这不起作用,以及为什么它的行为与第一个示例不同。我可以有你的建议来解决它吗?

标签: julia

解决方案


您需要Res就地更新。你现在做

Res = ones(S+S)

阴影输入Res。直接更新就好了Res


推荐阅读