kivantium活動日記

プログラムを使っていろいろやります

Rustによるニュートン法の実装と可視化

Rustの勉強として、1変数関数の簡単なニュートン法を実装してgnuplotによる可視化を行いました。

可視化ライブラリ

GitHub - SiegeLord/RustGnuplot: A Rust library for drawing plots, powered by Gnuplot.

[Rust plot]で検索して一番上に出てきたこれを使いました。
正直gnuplotの表示は嫌いなのでmatplotlibやggplotのラッパーが欲しいです。

実装

cargo new newton --bin

で新しいプロジェクトを作り、Cargo.toml

[dependencies]
gnuplot = "0.0.24"

を追記しました。

 f(x):=x^2-2=0の解を求めるプログラムは以下の通りです。

extern crate gnuplot;
use gnuplot::*;

// 解きたい関数 f
fn f(x: f64) -> f64 {
    x*x - 2.0
}

// fの導関数
fn df(x: f64) -> f64 {
    2.0*x
}

fn main() {
    // 関数のプロットに使うx座標たち
    let mut t = vec![];
    for i in 0..100i32 {
        t.push(i as f64 * 0.1 - 5.0);
    }
    let mut fg = Figure::new();
    // axes2d()はmutable borrowを作るので後でshow()するには別スコープを作る必要がある
    {
        let axes = fg.axes2d();
        // x軸を表示
        axes.set_x_axis(true, &[]);
        // 表示範囲の指定
        axes.set_x_range(Fix(0.0), Fix(5.0));
        axes.set_y_range(Fix(-5.0), Fix(30.0));
        // fのプロット
        axes.lines(t.iter(), t.iter().map(|&x| f(x)), &[Color("red")]);

        // ニュートン法
        let mut x = 4.0; // 初期値
        while f(x).abs() > 1e-10 { // f(x)が十分小さくなるまで続ける
            // 関数値のプロット
            axes.points(&[x], &[f(x)], &[Color("blue"), PointSymbol('O')]); 
            // 関数値に縦線を引く
            axes.lines(&[x, x], &[-5.0, 30.0], &[Color("blue")]);
            // 接線のプロット
            axes.lines(t.iter(), t.iter().map(|&p| df(x)*(p-x)+f(x)), &[Color("black")]);
            // 値の更新
            x = x - f(x)/df(x);
        }
        println!("solution: {}", x);
    }
    fg.show();
}

cargo runの結果は

solution: 1.4142135623730951

で、以下のようなプロットが表示されました。
f:id:kivantium:20180210133822p:plain

Rustは初めて書くので変なところがあったら教えてください。

そのうち自動微分とか精度保証とかやっていきたいところ。

特定商取引法に定められた事項は請求により遅滞なく提供する