LoginSignup
10
6

More than 5 years have passed since last update.

Rustで数値計算:Newton法

Posted at

Rustで数値計算を行う練習として、とりあえず一次元のNewton法を実装してみよう:

f(x) = 0

を解くために次のイテレーションを行う:

x_{n+1} = x_n - \frac{f(x_n)}{f^\prime(x_n)}

これを実装するには関数fを与える必要がある。
関数を関数にわたす部分はClosuresの章に書いてある。

fn newton<T1, T2>(f: T1, df: T2, mut x: f64, eps: f64) -> f64
    where T1: Fn(f64) -> f64,
          T2: Fn(f64) -> f64
{
    let mut y = f(x);
    while y.abs() > eps {
        x -= y / df(x);
        y = f(x);
    }
    x
}

2つのclosureは別の型になるのでT1,T2は別のキーワードにする必要がある。
関数のポインタでなく、コンパイル時に型が決まっている。

とりえあず微分は自分で解析的にやるとして使ってみよう:

fn main() {
    let mut x = -5.0;
    let f = |x| x * (x - 1.0);
    let df = |x| 2.0*x - 1.0;
    let eps = 1e-6;
    x = newton(&f, &df, x, eps);
    println!("Result = {:e}", x);
}
Result = -7.324879859764758e-11

最近ASTを解析して自動的に微分演算を実装する自動微分が流行ってるらしいが(要出展)、
そんな高度な事できないので、数値微分してみよう:

fn main() {
    let mut x = -5.0;
    let dx = 1e-7;
    let f = |x| x * (x - 1.0);
    let df = |x| (f(x+dx)-f(x))/dx;
    let eps = 1e-6;
    x = newton(&f, &df, x, eps);
    println!("Result = {:e}", x);
}
Result = -7.238726858155113e-11

課題

let df = |x| (f(x+dx)-f(x))/dx;

の部分を

let df = num_diff(&f, dx);

のように書きたいのだが、num_diffの定義の方法がわからなかった。
ラムダ式では型を書かなくて良いので定義できたが、別個の関数にする場合の型がわからない。
Haskellっぽく簡単に書けないのかな?

num_diff :: (Float -> Float) -> Float -> (Float -> Float)

Closuresには戻り値に使う場合はBox使えって書いているので、std::rc::Rc使えば

use std::rc::Rc;
fn num_diff(f: Rc<Fn(f64) -> f64>, dx: f64) -> Rc<Fn(f64) -> f64> {
    Rc::new(move |x| (f(x + dx) - f(x)) / dx)
}

fn main() {
    let mut x = -5.0;
    let f = Rc::new(|x| x * (x - 1.0));
    let df = num_diff(f.clone(), 1e-7);
    let eps = 1e-6;
    x = newton(&*f, &*df, x, eps);
    println!("Result = {:e}", x);
}

って書けるけど、これでいいのかな。。。。。

10
6
3

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
10
6