Julia で簡単な Newton 法の実装

関数 f_1(x) = x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3] についての Newton 法を実装してみると以下のようになり短いコードで書ける。

面白い点としては、

1. C や Matlab なら 2*x[1]*x[2] と書かなければいけないコードが 2x[1]x[2] と * を省略して書ける
2. UTF-8 でファイルを作っていれば、∇も普通の文字として扱える (全角空白が紛れ込むと、かなり厄介だけど)

実行するには
julia> include("newton.jl")
とする。
ちなみに、このコードなら短いので 15 分程度で書ける。

==== newton.jl ====

# -*- coding:utf-8 mode:Julia -*-

function f_1(x)
    return x[1]^4 - 2x[1]x[2]+ 2x[2]^2 - 4x[1] + 2x[3]
end
function ∇f_1(x)
  return [4x[1]^3-2x[2]-4; -2x[1]+4x[2]+2]
#    return [4*x[1]^3-2*x[2]-4; -2*x[1]+4*x[2]+2]
end
function ∇2f_1(x)
    return [12x[1]^2 -2; -2 4]
end


ε = 1.0e-8
k = 0
x = [1; 1]
norm_∇f = 0.0
while k < 100
    ∇f = ∇f_1(x)
    norm_∇f = norm(∇f)
    if norm_∇f < ε
        println("Convergence!!")
        break
    end
    ∇2f = ∇2f_1(x)
    @printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f
    x = x - ∇2f \ ∇f
    k = k+1
end
@printf "k = %d : x = [%.2e, %.2e], norm(∇f) = %.2e\n" k x[1] x[2] norm_∇f

==========