３年振りのRustなので数値計算で構文を復習してみた

この記事は Rustその2 Advent Calendar 2019 5日目の記事です。

はじめに

こんにちは、こんばんは。kaneshin です。Rustも2018 Editionとなっていますが、自分が最後にRustを触ったのが3年前の2016年となっていて、基本構文さえも覚束ないので初心者の人がRustらしくコードを変更させていく記事にさせてもらいます。（本当はネタが思いつかなかった）

普段は Go をよく書いているので、言語仕様が豊富な言語はワクワクしますね！

数値積分のコードを最適化する

さて、題材として今回は数値積分を取り上げました。

$\displaystyle{ \int_{-2}^4 (x^2-x-1) dx = 12 }$

www.wolframalpha.com

この積分式を数値積分の解法である Midpoint Rule （中点法）と Trapezoid Rule （台形公式）の二つを利用します。

台形公式 - Wikipedia

ひとまず何も考えずにコードを書く

すべて32bitで数値解析を行います。

fn f(x: f32) -> f32 {
    return x * x - x - 1.0;
}

fn midpoint_rule(a: f32, b: f32, n: i32) -> f32 {
    let h = (b - a) / n as f32;
    let mut r = 0.0;
    let mut x = a + h / 2.0;

    for _ in 0..n {
        r += f(x);
        x += h;
    }
    return r * h;
}

fn trapezoid_rule(a: f32, b: f32, n: i32) -> f32 {
    let h = (b - a) / n as f32;
    let mut r = 0.0;
    let mut x = a + h;
    let mut y1 = f(a);

    for _ in 0..n {
        let y2 = f(x);
        r += y1 + y2;
        x += h;
        y1 = y2;
    }
    return r * h / 2.0;
}

fn main() {
    println!("Midpoint Rule: {}", midpoint_rule(-2.0, 4.0, 1000));
    println!("Trapezoid Rule: {}", trapezoid_rule(-2.0, 4.0, 1000));
}

これの実行結果は下記のように出力されます。

Midpoint Rule: 12.000174
Trapezoid Rule: 12.000226

この関数だと中点法の方が精度高そうですね。

いろいろ施してみた

やったこととしては：

trait と impl を使って抽象化させた
- 関数で数値解析をしていたが、 impl でそれぞれの構造体に関数を定義
- 結果出力は trait で定義している
- fn foo(integral: impl Integral); のようにして呼び出せる
Tuple を使って積分区間を表現した

trait Integral {
    fn method(&self) -> String;
    fn calc(&self) -> f32;
    fn print(&self) {
        println!("{}: {}", self.method(), self.calc());
    }
}

struct Interval(f32, f32);

struct MidpointRule {
    function: fn(f32) -> f32,
    interval: Interval,
    step: i32,
}

impl Integral for MidpointRule {
    fn method(&self) -> String {
        return String::from("Midpoint Rule");
    }

    fn calc(&self) -> f32 {
        let a = self.interval.0;
        let b = self.interval.1;
        let n = self.step;
        let h = (b - a) / n as f32;
        let mut r = 0.0;
        let mut x = a + h / 2.0;

        for _ in 0..n {
            r += (self.function)(x);
            x += h;
        }
        return r * h;
    }
}

struct TrapezoidRule {
    function: fn(f32) -> f32,
    interval: Interval,
    step: i32,
}

impl Integral for TrapezoidRule {
    fn method(&self) -> String {
        return String::from("Trapezoid Rule");
    }

    fn calc(&self) -> f32 {
        let a = self.interval.0;
        let b = self.interval.1;
        let n = self.step;
        let h = (b - a) / n as f32;
        let mut r = 0.0;
        let mut x = a + h;
        let mut y1 = (self.function)(a);

        for _ in 0..n {
            let y2 = (self.function)(x);
            r += y1 + y2;
            x += h;
            y1 = y2;
        }
        return r * h / 2.0;
    }
}

fn main() {
    fn f(x: f32) -> f32 {
        return x * x - x - 1.0;
    }

    MidpointRule {
        function: f,
        interval: Interval(-2.0f32, 4.0f32),
        step: 1000,
    }.print();

    TrapezoidRule {
        function: f,
        interval: Interval(-2.0f32, 4.0f32),
        step: 1000,
    }.print();
}