
using LinearAlgebra,IntervalArithmetic, ForwardDiff
X = [0.707107 ± 1e-8, 0.707107 ± 1e-8]
c = mid.(X)
f(x, y) = x^2 + y^2 - 1
g(x, y) = x - y
F( (x, y) ) = [f(x, y); g(x, y)]
DF = ForwardDiff.jacobian(F,c)
R = inv(DF)
s = X - R*F(X)
@show s
@show s .⊂ X;

s = Interval{Float64}[[0.707106, 0.707107], [0.707106, 0.707107]]
s .⊂ X = Bool[0, 0]


using LinearAlgebra,IntervalArithmetic, ForwardDiff
X = [(0.8..0.9),(-1..(-0.5))]
f(x, y) = x^2 + y^2 - 1
g(x, y) = x^3 + y^4
h( (x, y) ) = [f(x, y); g(x, y)]
# ForwardDiff.jacobian(g, X::IntervalBox) = ForwardDiff.jacobian(g, X.v)
J = ForwardDiff.jacobian(h,X)

2×2 Matrix{Interval{Float64}}:
 [1.59999, 1.80001]  [-2, -1]
 [1.91999, 2.43001]   [-4, -0.5]


h(X)

2-element Vector{Interval{Float64}}:
 [-0.110001, 0.810001]
  [0.574499, 1.72901]


using LinearAlgebra,IntervalArithmetic, ForwardDiff
# 区間Xを定義
X = [(0.6..0.8),(0.6..0.8)]
# 計算対象となる方程式を定義
f(x, y) = x^2 + y^2 - 1
g(x, y) = x - y
F( (x, y) ) = [f(x, y); g(x, y)]
# ForwardDiff.jacobian(g, X::IntervalBox) = ForwardDiff.jacobian(g, X.v)
# 区間Xにおけるヤコビ行列を計算
iDF = ForwardDiff.jacobian(F,X)

2×2 Matrix{Interval{Float64}}:
     [1.19999, 1.60001]        [1.19999, 1.60001]
 [1, 1]                  [-1, -1]


c = mid.(X)
ic = map(Interval,c)
# Rを計算するためのf'(c)を求める
DF = ForwardDiff.jacobian(F,c) # 区間演算なし

2×2 Matrix{Float64}:
 1.4   1.4
 1.0  -1.0


# 区間Xにおけるヤコビ行列Jの逆行列Rを定義する
R = inv(DF) # 区間演算なし

2×2 Matrix{Float64}:
 0.357143   0.5
 0.357143  -0.5


M = Matrix{Float64}(I,size(R)) - R*iDF

2×2 Matrix{Interval{Float64}}:
 [-0.0714286, 0.0714286]  [-0.0714286, 0.0714286]
 [-0.0714286, 0.0714286]  [-0.0714286, 0.0714286]


R*F(ic)

2-element Vector{Interval{Float64}}:
 [-0.00714286, -0.00714285]
 [-0.00714286, -0.00714285]


# Krawczyk写像を計算
K = c - R*F(ic) + M*(X - c)　# 区間演算必要

2-element Vector{Interval{Float64}}:
 [0.692857, 0.721429]
 [0.692857, 0.721429]


# 収束判定
K .⊂ X

2-element BitVector:
 1
 1


function krawczyk(F,X)
    iDF = ForwardDiff.jacobian(F,X)
    c = mid.(X); ic = map(Interval,c)
    DF = ForwardDiff.jacobian(F,c)
    R = inv(DF)
    M = Matrix{Float64}(I,size(R)) - R*iDF
    return c - R*F(ic) + M*(X - c)
end

krawczyk (generic function with 1 method)


K = krawczyk(F,X)
tol = 5e-10
while maximum(radius,K) >= tol
    K = krawczyk(F,K)
   @show radius.(K)
end
K

radius.(K) = [0.00028860028860044906, 0.00028860028860044906]
radius.(K) = [1.1779002662137827e-7, 1.1779002662137827e-7]
radius.(K) = [1.9761969838327786e-14, 1.9761969838327786e-14]

2-element Vector{Interval{Float64}}:
 [0.707106, 0.707107]
 [0.707106, 0.707107]


err = 0.7;
X = [0.7 ± err , 0.7 ± err ];
c = mid.(X);
K = krawczyk(F,X);
K .⊂ X

2-element BitVector:
 0
 0


# rを計算
r = abs.(R*F(c))

2-element Vector{Float64}:
 0.00714285714285707
 0.00714285714285707


# uを計算
u = r .+ (sum(r)/length(r))

2-element Vector{Float64}:
 0.01428571428571414
 0.01428571428571414


# 候補区間Xを新たに定める
X_new = c .± u;
X = X_new;
K = krawczyk(F,X);
K .⊂ X

2-element BitVector:
 1
 1


using LinearAlgebra, IntervalArithmetic, ForwardDiff

# 候補区間Xを定義
X = [(0.6..0.7),(0.6..0.8)]

# 計算対象となる方程式を定義
f(x, y) = x^2 + y^2 - 1
g(x, y) = x^2 - y^4
F( (x, y) ) = [f(x, y); g(x, y)]

F (generic function with 1 method)


#ニュートン法で近似解を計算する
function newton(F,x0)
    #初期値を設定
    tol = 5e-10; count = 0;
    x = x0;
    Fx = F(x);
    #条件の範囲内で計算を回す
    while maximum(abs,Fx) ≥ tol && count ≤ 20
        DF = ForwardDiff.jacobian(F,x);
        x -= DF\Fx;
        Fx = F(x);
        count += 1;
    end
    return x
end

#クラフチック法を計算する
function krawczyk(F,X)
    iDF = ForwardDiff.jacobian(F,X);
    c = mid.(X); ic = map(Interval,c);
    DF = ForwardDiff.jacobian(F,c);
    R = inv(DF);
    M = Matrix{Float64}(I,size(R)) - R*iDF;
    #クラフチック写像の値を返す
    return c - R*F(ic) + M*(X - c)
end

#最終的に完成した関数
function verifynlss_krawczyk(F,c)
    DF = ForwardDiff.jacobian(F,c)
    R = inv(DF)
    r = abs.(R*F(c))
    u = r .+ (sum(r)/length(r))
    X = c .± u
    K = krawczyk(F,X)
    #範囲内に入っていたら、さらに範囲の精度をあげていく
    if all(K .⊂ X)
        tol = 5e-10
        count = 0
        while maximum(radius,K) >= tol && count ≤ 100
            K = krawczyk(F,K)
            count += 1
        end
        success = 1
        return success, K
    end
    println("Oh my way, verification is failed...return a improved approximate solution") # cをNewton法で改善しても良い。
    success = 0
    return success, newton(F,c)
end

success, X = verifynlss_krawczyk(F,[0.61,0.78])
if success == 0
    success, X = verifynlss_krawczyk(F,X)
end
success, X

(1, Interval{Float64}[[0.618033, 0.618034], [0.786151, 0.786152]])


# include("IntervalLinearAlgebra.jl");

#区間連立1次方程式を解く関数
function verifylss_iAib(iA,ib) 
    A = mid.(iA)
    b = mid.(ib)
    x̄ = A\b
    n = length(x̄)
    R = inv(A)
    #########
    G = Matrix{Float64}(I, n, n) - R*iA
    α = opnorm(G,Inf)# Interval arithmetic
    #########
    if α < 1
        x̄ = map(Interval,x̄)
        r = iA*x̄ - ib # Interval arithmetic
        Rr = R*r
        err = abs.(Rr) + supremum(norm(Rr,Inf))/(1-α)*(abs.(G)*ones(n)) # Interval arithmetic
    else
        println("Oh my way, verification is failed...")
        err = nan
    end
    return x̄ .± supremum.(err)
end

verifylss_iAib (generic function with 1 method)


using LinearAlgebra,IntervalArithmetic, ForwardDiff
# 区間Xを定義
X = [(0.6.. 0.7),(0.7.. 0.8)]
c = mid.(X)
ic = map(Interval,c)

# 計算対象となる方程式を定義
f(x, y) = x^2 + y^2 - 1
g(x, y) = x^2 - y^4
F( (x, y) ) = [f(x, y); g(x, y)]

# 区間Xにおけるヤコビ行列を計算
iM  = ForwardDiff.jacobian(F,X) # 区間演算必要
ifc = F(ic)

2-element Vector{Interval{Float64}}:
 [-0.0150001, -0.0149999]
  [0.106093, 0.106094]


#区間連立1次方程式を計算
verifylss_iAib(iM,ifc)

2-element Vector{Interval{Float64}}:
  [0.0206957, 0.0432826]
 [-0.047109, -0.0283388]


#N(c,X)を計算
N = ic - verifylss_iAib(iM,ifc)

2-element Vector{Interval{Float64}}:
 [0.606717, 0.629305]
 [0.778338, 0.797109]


#収束判定
N .⊂ X

2-element BitVector:
 1
 1


#ニュートン法を計算
function newton(F,x0)
    #初期値を設定
    tol = 5e-10; count = 0;
    x = x0;
    Fx = F(x);
    #条件によってニュートン法をまわす
    while maximum(abs,Fx) ≥ tol && count ≤ 20
        DF = ForwardDiff.jacobian(F,x);
        x -= DF\Fx;
        Fx = F(x);
        count += 1;
    end
    return x
end

#N(c,X)を計算する関数
function IntervalNewton(F,X)
    c = mid.(X);
    ic = map(Interval,c);
    iM = ForwardDiff.jacobian(F,X);
    ib = F(ic);
    #N(c,X)の値を返す
    return ic - verifylss_iAib(iM,ib)
#     return ic - iM\ib
end

#最終的に構築した関数
function verifynlss_IntervalNewton(F, c)
    DF = ForwardDiff.jacobian(F,c);
    R = inv(DF);
    r = abs.(R*F(c));
    u = r .+ (sum(r)/length(r));
    X = c .± u;
    K = IntervalNewton(F,X);
    #範囲内に入っていたら、さらに解の精度をあげる
    if all(K .⊂ X)
        tol = 5e-10;
        while maximum(radius,K) >= tol
            K = IntervalNewton(F,K)
        end
        success = 1
        return success, K
    end
    println("Oh my way, verification is failed...return a improved approximate solution") 
    success = 0
    return success, newton(F,c)
end

verifynlss_IntervalNewton (generic function with 1 method)


verifynlss_IntervalNewton(F,[0.7,0.8])

(1, Interval{Float64}[[0.618033, 0.618034], [0.786151, 0.786152]])


using LinearAlgebra,IntervalArithmetic, ForwardDiff

#\lambdaを設定
lam = 3.82843

#解く方程式を設定
f(x, y, z) = x - lam*z*(1-z) 
g(x, y, z) = y - lam*x*(1-x)
h(x, y, z) = z - lam*y*(1-y)
F( (x, y, z) ) = [f(x, y, z); g(x, y, z); h(x, y, z)]

#候補区間を設定
X = [(0.9.. 1.0),(0.1.. 0.2),(0.5.. 0.6)]

3-element Vector{Interval{Float64}}:
 [0.899999, 1]
 [0.0999999, 0.200001]
 [0.5, 0.600001]


x0 = newton(F,mid.(X)) # 区間演算なし

3-element Vector{Float64}:
 0.9562724713863567
 0.16008745377675246
 0.5147686339721098


using BenchmarkTools
@btime verifynlss_IntervalNewton($F,$x0)

  30.000 μs (142 allocations: 12.19 KiB)

(1, Interval{Float64}[[0.956272, 0.956273], [0.160087, 0.160088], [0.514768, 0.514769]])


@btime verifynlss_krawczyk($F,$x0)

  23.541 μs (144 allocations: 12.14 KiB)

(1, Interval{Float64}[[0.956272, 0.956273], [0.160087, 0.160088], [0.514768, 0.514769]])

Juliaで非線型方程式の解の精度保証付き数値計算¶

Krawczyk（クラフチック）法¶

縮小写像の原理¶

Krawczyk写像と簡易ニュートン写像の関係¶

平均値形式¶

自動微分を利用したヤコビ行列の計算¶

Krawczyk法の実装¶

区間ニュートン法¶

Krawczyk法との違い¶

区間ニュートン法の実装¶

参考文献¶