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);
}
って書けるけど、これでいいのかな。。。。。