3.2 Julia でNewton法

まず Julia をインストールしておこう。 「インストール」

newton.jl

# newton.jl

using Printf
using Plots

function newton(f, df, x0, eps)
  x = x0
  for i=1:10
    dx = f(x) / df(x)
    x = x - dx
    #@printf("Δx=%e, x=%g, f(x)=%e\n", dx, x, f(x))
    if abs(dx) < eps
      return x
    end
  end
  println("newton: 収束しませんでした。修正量=$(dx)")
  return x
end

#
println("問題1")
f(x) = x^3 - 2x - 5
df(x)=3x^2-2

println("通常の浮動小数点演算")
println(newton(f, df, 2.0, 1e-14))

println("BigFloat を用いた10進100桁計算")
# 2進350桁(10進100桁のためには330桁あれば良い)
setprecision(350)
println(newton(f, df, BigFloat(2), 1e-100))

#
println("問題2")
f(x) = tan(x)+x
df(x)=1/cos(x)^2+1
println("通常の浮動小数点演算")
println(newton(f, df, 2.0, 1e-14))

println("10進100桁計算")
println(newton(f, df, BigFloat(2), 1e-100))

@printf() の注釈を外すと中間結果が見える。

Julia でNewton法を試す (% や julia> の右側を入力する)

% curl -O http://nalab.mind.meiji.ac.jp/~mk/misc/20211210/newton.jl

% julia
julia> include("newton.jl")
通常の浮動小数点演算
2.0945514815423265
BigFloat を用いた10進100桁計算
2.0945514815423265914823865405793029638573061056282391803041285290453121899834836671462672817771577578608389
通常の浮動小数点演算
2.028757838110434
10進100桁計算
2.0287578381104342235769711247347143761083800287593940888171660744498665031042762345922795150425630639023986

あるいは

% julia newton.jl

このプログラムを発展させて、 解 $u(x,t)$ を計算してグラフを描くプログラムを描いてみよう。