↑Programmieren in Rust

Algorithmen: Numerik

Inhaltsverzeichnis

1. Bisektionsverfahren
2. Numerische Differentiation
3. Numerische Integration
4. Interpolation

Bisektionsverfahren

fn bisection(
    mut a: f64, mut b: f64, epsilon: f64,
    f: &dyn Fn(f64) -> f64
) -> f64 {
    for _ in 0..100 {
        let c = 0.5*(a + b);
        if f(c) == 0.0 {return c;}
        if f(a).signum() != f(c).signum()  {b = c;} else {a = c;}
        if (b - a).abs() < epsilon {break;}
    }
    0.5*(a + b)
}

fn checked_bisection(
    mut a: f64, mut b: f64, epsilon: f64, nmax: u32,
    f: &dyn Fn(f64) -> f64
) -> Option<f64> {
    if !(f(a) < 0.0 && f(b) > 0.0 || f(a) > 0.0 && f(b) < 0.0) {
        return None;
    }
    for _ in 0..nmax {
        let c = 0.5*(a + b);
        if f(c) == 0.0 {return Some(c);}
        if f(a).signum() != f(c).signum()  {b = c;} else {a = c;}
        if (b - a).abs() < epsilon {return Some(0.5*(a + b));}
    }
    None
}

Numerische Differentiation

Erste Ableitung

fn diff(f: &dyn Fn(f64) -> f64, x: f64) -> f64 {
    let h = 0.001;
    (f(x+h) - f(x-h))/(2.0*h)
}

Zweite Ableitung

fn diff2(f: &dyn Fn(f64) -> f64, x: f64) -> f64 {
    let h = 0.001;
    (f(x+h) - 2.0*f(x) + f(x-h))/(h*h)
}

Dritte Ableitung

fn diff3(f: &dyn Fn(f64) -> f64, x: f64) -> f64 {
    let h = 0.004;
    (0.5*f(x+2.0*h) - f(x+h) + f(x-h) - 0.5*f(x-2.0*h))/(h*h*h)
}

Numerische Integration

Simpsonregel

fn simpson(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: u32) -> f64 {
    let h = (b - a)/(n as f64);
    let mut s = 0.0;
    for i in 0..n {
        let x = a + h*(i as f64);
        s += f(x) + 4.0*f(x + 0.5*h) + f(x + h);
    }
    s*h/6.0
}

Gauß-Legendre-Quadratur

static GL32: [[f64; 2]; 32] = [
    [-0.9972638618494816, 0.007018610009470141],
    [-0.9856115115452684, 0.01627439473090563],
    [-0.9647622555875064, 0.02539206530926208],
    [-0.9349060759377397, 0.03427386291302148],
    [-0.8963211557660522, 0.04283589802222668],
    [-0.8493676137325699, 0.05099805926237622],
    [-0.7944837959679425, 0.05868409347853551],
    [-0.7321821187402897, 0.06582222277636184],
    [-0.6630442669302152, 0.07234579410884850],
    [-0.5877157572407623, 0.07819389578707028],
    [-0.5068999089322295, 0.08331192422694673],
    [-0.4213512761306354, 0.08765209300440380],
    [-0.3318686022821277, 0.09117387869576390],
    [-0.2392873622521371, 0.09384439908080457],
    [-0.1444719615827965, 0.09563872007927489],
    [-0.04830766568773831,0.09654008851472778],
    [ 0.04830766568773831,0.09654008851472778],
    [ 0.1444719615827965, 0.09563872007927489],
    [ 0.2392873622521371, 0.09384439908080457],
    [ 0.3318686022821277, 0.09117387869576390],
    [ 0.4213512761306354, 0.08765209300440380],
    [ 0.5068999089322295, 0.08331192422694673],
    [ 0.5877157572407623, 0.07819389578707028],
    [ 0.6630442669302152, 0.07234579410884850],
    [ 0.7321821187402897, 0.06582222277636184],
    [ 0.7944837959679425, 0.05868409347853551],
    [ 0.8493676137325699, 0.05099805926237622],
    [ 0.8963211557660522, 0.04283589802222668],
    [ 0.9349060759377397, 0.03427386291302148],
    [ 0.9647622555875064, 0.02539206530926208],
    [ 0.9856115115452684, 0.01627439473090563],
    [ 0.9972638618494816, 0.007018610009470141]
];

fn gauss(f: &dyn Fn(f64) -> f64, a: f64, b: f64, n: u32) -> f64 {
    let h = (b - a)/(n as f64);
    let p = 0.5*h;
    let mut s = 0.0;
    for j in 0..n {
        let q = p + a + (j as f64)*h;
        let mut sj = 0.0;
        for t in &GL32 {
            sj += t[1]*f(p*t[0] + q);
        }
        s += sj;
    }
    p*s
}

Interpolation

Stückweise lineare interpolation

// Für äquidistante Stützstellen [x0 + k*d, y[k]].

fn pli(x0: f64, d: f64, y: Vec<f64>) -> impl Fn(f64) -> f64 {
    let n = y.len() as isize;
    move |x| {
        let k = f64::floor((x-x0)/d) as isize;
        if k < 0 || k+1 >= n {f64::NAN}
        else {
            let k = k as usize;
            y[k] + (y[k+1] - y[k])/d*(x - x0 - (k as f64)*d)
        }
    }
}