32 Julia で遊ぶ (1)

3年くらい前に一度 Julia で遊んだことがあったのだけど、 しばらく忘れていて、 先日某先生に「Julia はどうですか」みたいなことを聞かれて、 今日また別件で目にしたので、久しぶりにネットを見に行く。


SIAM Review (有名なジャーナルです) で紹介されたりしている?ほー。

J. Bezanson, A. Edelman, S. Karpinski, V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review, Vol. 59, pp. 65-98 (2017).
(これは学内ネットワークからは読める。)

Juliaのインストール     Homebrew では
brew update
brew cask install julia
とせよ、と言われている。

MacPorts では
sudo port install julia
かな。

でも、今回は、 Julialang から入手する。 今だと julia-1.2.0-mac64.dmg を保存して、ダブルクリックすると

図 1: Julia-1.2 を Applications までドラッグして入れればOK
Image Julia-install

これで、 アプリケーション ディレクトリィに、 Julia-1.2.app と言うディレクトリィが出来る。ダブル・クリックすると、 ターミナルが起動してその中で Julia が起動する。 CUI なのか。では、ターミナルや iTerm2 から呼び出すとして、 エイリアスの登録でもしておくか。
alias julia=/Applications/Julia-1.2.app/Contents/Resources/julia/bin/julia
とするのかな?

$ julia
               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.2.0 (2019-08-20)
 _/ |\__'_|_|_|\__'_|  |  Official https://julialang.org/ release
|__/                   |

julia>

今回は「チュートリアル」で紹介されていたオープン・テキストの Giray Ökten, First Semester in Numerical Analysis with Julia, 2018 をパラパラしながら試す。

ベクトル     ベクトルの文法は MATLAB ライクだけれど独特なところもある。
julia> x=[1,2,3,4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> y=[1;2;3;4]
4-element Array{Int64,1}:
 1
 2
 3
 4

julia> z=[1.0;2;3;4]
4-element Array{Float64,1}:
 1.0
 2.0
 3.0
 4.0
どれも縦ベクトル。,; で区切る。混ぜられない。 すべて Int64 か、すべて Float64 か。 Int32 ではないのか。まー、確かに 64 ビットの時代だよね。

ふと、多倍長はどうなるのかなあ?と気になった。
ちょっと脱線
GNU MPFR が入っていて、任意精度演算と言うのが出来るんだ。
julia> setprecision(300)
300
julia> x=BigFloat(1)
1.0
julia> x/3
0.33333333333333333333333333333333333333333333333333333333333333333333333333333333333333333342
setprecision() はビット数の指定をするんだ。 そのうちスピードを調べたい。

(調子に乗って) 区間演算はどうだろう? MATLAB は INTLAB のようなツールボックスを可能にする言語仕様だったのだけれど、 Julia はどうか? …すでに色々チャレンジしている人がいるんだ (https://swim2019.ensta-paris.fr/papers/vorontsova.pdf)。へー。 今度試してみよう。

MATLABのように ' で転置(正確には Hermite 共役)される。
julia> x'*x
30

こういうこともできる。
julia> x=[10*i for i=1:5]
5-element Array{Int64,1}:
 10
 20
 30
 40
 50
Python調?

julia> first(x)
10
julia> last(x)
50
julia> minimum(x)
10
julia> maximum(x)
50
julia> sum(x)
150
julia> append!(x,60)
6-element Array{Int64,1}:
 10
 20
 30
 40
 50
 60
julia> sin.(x)
6-element Array{Float64,1}:
 -0.5440211108893698
  0.9129452507276277
 -0.9880316240928618
  0.7451131604793488
 -0.26237485370392877
 -0.3048106211022167

グラフ描き    

julia> と言うプロンプトに対して、] と打つと、 package mode になる。add PyPlot とすると、 ネットから必要なファイルをダウンロードして、 インストーする。結構時間がかかるね。 package mode から元に戻るには backspace (control+H) を打つ。

julia> using PyPlot
[ Info: : Precompiling PyPlot [d330b81b-6aea-500a-939a-2ce795aea3ee]
むにゃむにゃ
色々なものをダウンロードし始めた。

  cycler             pkgs/main/osx-64::cycler-0.10.0-py37_0
  freetype           pkgs/main/osx-64::freetype-2.9.1-hb4e5f40_0
  kiwisolver         pkgs/main/osx-64::kiwisolver-1.1.0-py37h0a44026_0
  libpng             pkgs/main/osx-64::libpng-1.6.37-ha441bb4_0
  matplotlib         pkgs/main/osx-64::matplotlib-3.1.1-py37h54f8f79_0
  pyparsing          pkgs/main/noarch::pyparsing-2.4.2-py_0
  python-dateutil    pkgs/main/osx-64::python-dateutil-2.8.0-py37_0
  pytz               pkgs/main/noarch::pytz-2019.3-py_0
  tornado            pkgs/main/osx-64::tornado-6.0.3-py37h1de35cc_0

あ、PyPlot って、Python の Plot ということか。

(Catalina では何も言われないが、High Sierra や El Capitan では
┌ Warning: PyPlot is using tkagg backend, which is known to cause crashes on MacOS (#410); use the MPLBACKEND environment variable to request a different backend.
└ @ PyPlot ~/.julia/packages/PyPlot/4wzW1/src/init.jl:192
という警告が表示された。Tkagg ではなくQt4Agg というやつか? ~/.matplotlib/matplotlibrcbackend : Qt4Agg と書くと警告は出なかった。 これまでは backend : TkAgg と書いてあった。)

julia> というプロンプトに対して
using PyPlot
x = range(0,stop=2*pi,length=1000)
y = sin.(3x);
plot(x, y, color="red", linewidth=2.0, linestyle="--")
title("The sine function");

Image julia_figure_1
フロッピー・ディスク (まだこれを使うのか) のアイコンをクリックすると、 保存できる。

(range() って何?
julia> range(0,stop=1,length=10)
0.0:0.1111111111111111:1.0

julia> range(0,stop=1,length=11)
0.0:0.1:1.0
$ [a,b]$$ n$ 等分するなら range(a,stop=b,length=n+1) かな。 MATLAB で linspace(a,b,n+1) とするのと同じだ。

ところでどうやって終了するのか? 当てずっぽうで quit, exit, bye を試したけれどダメ。 control+D で終了できたけれど、まさかねえ。 真面目に調べたら (「Julia練習帳」)、 exit() とか exit(1) とかするんだって。

プログラム

拡張子は .jl なんだ。
bisection.jl
function bisection(f::Function,a,b,eps,N)
    n=1
    p=0. # to ensure the value of p carries out of the while loop
    while n<=N
        p = a+(b-a)/2
        if f(p)==0 || abs(a-b)<eps
            return println("p is $p and the iteration number is $n")
        end
        if f(a)f(p)<0
            b=p
        else
            a=p
        end
        n=n+1
    end
    y=f(p)
    println("Method did not converge. The last iteration gives $p with function value $y")
end
ここは由緒正しい Bombelli の方程式 $ x^3-15x-4=0$ で試そう。
julia> include("bisection.jl")
bisection (generic function with 1 method)

julia> bisection(x->x^3-15x-4,2,6,10^(-15),20)
p is 4.0 and the iteration number is 1

julia> bisection(x->x^3-15x-4,2,5,10^(-15),20)
Method did not converge. The last iteration gives 4.000000953674316 with function value 3.147126335534267e-5

julia> bisection(x->x^3-15x-4,2,5,10^(-15),30)
Method did not converge. The last iteration gives 4.000000000931323 with function value 3.073364496231079e-8

julia> bisection(x->x^3-15x-4,2,5,10^(-15),100)
p is 4.0 and the iteration number is 52

julia> function g(x)
       return x^3-15x-4
       end
g (generic function with 1 method)

julia> bisection(g,2,5,10^(-15),100)
p is 4.0 and the iteration number is 52
おお、面白い (今さらながら二分法は遅いね)。

文字列の引用符は、クォーテーション・マーク " なんだ。


ゼミで、とりあえず C 言語の実習を始めたけれど、 いつまでも C べったりでもないだろう、 Python でもやろうかな?と思っていたけれど、 案外 Julia も面白いかもしれないな。


多分続く。


(2019/11/6) アマゾンで書籍を検索してみて、 アダルト扱いになってしまう (苦笑)。 ヒットした和書は1つだけ。

Juliaプログラミングクックブック ―言語仕様からデータ分析, 機械学習、数値計算まで
アマゾンではレビューなし。Amazon では星3つ。 読者を選ぶ、みたいな微妙な評価だったけれど、 どちらかというと、自分は想定されている読者のような気がしたので買うことに。

翻訳者が書いた文章もついているけれど、 その部分は微妙な感じがした。

全体的に人に勧める感じの本ではないのかな。

(2019/11/18) これで数値解析入門のレクチャーをするとどうなるか? と考えて、サンプル・プログラムを書き始めた。 面白いのは、二分法のプログラムを書いて、 ふと「引数に BigFloat() を渡したら?」と思いついて実行したら、 多倍長計算をしてくれたこと。なるほど。

(2019/11/17) パッケージを加えるのに、次のようなやり方を指示された。 ] 打ってパッケージモードで作業するのと違いはあるのか?
using Pkg
Pkg.add("DifferentialEquations")
Pkg.add("Plots")
using DifferentialEquations
using Plots

(2019/11/28) Julia で Plots のインストールの途中でひっかかると、後が面倒。
julia> using Plots
[ Info: Precompiling Plots [91a5bcdd-55d7-5caf-9e0b-520d859cae80]
ERROR: LoadError: UndefVarError: ffmpeg not defined
Stacktrace:
....
https://github.com/JuliaIO/FFMPEG.jl/issues/12 に書いてあるように、 ] としてから

(v1.2) pkg> add FFMPEG#master; build FFMPEG; test FFMPEG

桂田 祐史
2020-04-20