refactor newton's method code
This commit is contained in:
parent
cc2722e58c
commit
3e00657ade
124
src/function.rs
124
src/function.rs
@ -2,7 +2,7 @@
|
|||||||
|
|
||||||
use crate::function_output::FunctionOutput;
|
use crate::function_output::FunctionOutput;
|
||||||
#[allow(unused_imports)]
|
#[allow(unused_imports)]
|
||||||
use crate::misc::{debug_log, SteppedVector};
|
use crate::misc::{debug_log, newtons_method, SteppedVector};
|
||||||
|
|
||||||
use crate::egui_app::{DEFAULT_FUNCION, DEFAULT_RIEMANN};
|
use crate::egui_app::{DEFAULT_FUNCION, DEFAULT_RIEMANN};
|
||||||
use crate::parsing::BackingFunction;
|
use crate::parsing::BackingFunction;
|
||||||
@ -228,109 +228,35 @@ impl FunctionEntry {
|
|||||||
// Finds roots
|
// Finds roots
|
||||||
fn roots(&mut self) {
|
fn roots(&mut self) {
|
||||||
let resolution: f64 = (self.pixel_width as f64 / (self.max_x - self.min_x).abs()) as f64;
|
let resolution: f64 = (self.pixel_width as f64 / (self.max_x - self.min_x).abs()) as f64;
|
||||||
let mut root_list: Vec<Value> = Vec::new();
|
self.output.roots = Some(
|
||||||
let mut last_ele: Option<Value> = None;
|
newtons_method(
|
||||||
for ele in self.output.back.as_ref().unwrap().iter() {
|
resolution,
|
||||||
if last_ele.is_none() {
|
self.min_x..self.max_x,
|
||||||
last_ele = Some(*ele);
|
self.output.back.to_owned().unwrap(),
|
||||||
continue;
|
&|x: f64| self.function.get(x),
|
||||||
}
|
&|x: f64| self.function.get_derivative_1(x),
|
||||||
|
)
|
||||||
let last_ele_signum = last_ele.unwrap().y.signum();
|
.iter()
|
||||||
let ele_signum = ele.y.signum();
|
.map(|x| Value::new(*x, self.function.get(*x)))
|
||||||
|
.collect(),
|
||||||
if last_ele_signum.is_nan() | ele_signum.is_nan() {
|
);
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
if last_ele_signum != ele_signum {
|
|
||||||
// Do 50 iterations of newton's method, should be more than accurate
|
|
||||||
let x = {
|
|
||||||
let mut x1: f64 = last_ele.unwrap().x;
|
|
||||||
let mut x2: f64;
|
|
||||||
let mut fail: bool = false;
|
|
||||||
loop {
|
|
||||||
x2 = x1 - (self.function.get(x1) / self.function.get_derivative_1(x1));
|
|
||||||
if !(self.min_x..self.max_x).contains(&x2) {
|
|
||||||
fail = true;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (x2 - x1).abs() < resolution {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
x1 = x2;
|
|
||||||
}
|
|
||||||
|
|
||||||
match fail {
|
|
||||||
true => f64::NAN,
|
|
||||||
false => x1,
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
if !x.is_nan() {
|
|
||||||
root_list.push(Value::new(x, self.function.get(x)));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
last_ele = Some(*ele);
|
|
||||||
}
|
|
||||||
self.output.roots = Some(root_list);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
// Finds extrema
|
// Finds extrema
|
||||||
fn extrema(&mut self) {
|
fn extrema(&mut self) {
|
||||||
let resolution: f64 = (self.pixel_width as f64 / (self.max_x - self.min_x).abs()) as f64;
|
let resolution: f64 = (self.pixel_width as f64 / (self.max_x - self.min_x).abs()) as f64;
|
||||||
let mut extrama_list: Vec<Value> = Vec::new();
|
self.output.extrema = Some(
|
||||||
let mut last_ele: Option<Value> = None;
|
newtons_method(
|
||||||
for ele in self.output.derivative.as_ref().unwrap().iter() {
|
resolution,
|
||||||
if last_ele.is_none() {
|
self.min_x..self.max_x,
|
||||||
last_ele = Some(*ele);
|
self.output.derivative.to_owned().unwrap(),
|
||||||
continue;
|
&|x: f64| self.function.get_derivative_1(x),
|
||||||
}
|
&|x: f64| self.function.get_derivative_2(x),
|
||||||
|
)
|
||||||
let last_ele_signum = last_ele.unwrap().y.signum();
|
.iter()
|
||||||
let ele_signum = ele.y.signum();
|
.map(|x| Value::new(*x, self.function.get(*x)))
|
||||||
|
.collect(),
|
||||||
if last_ele_signum.is_nan() | ele_signum.is_nan() {
|
);
|
||||||
continue;
|
|
||||||
}
|
|
||||||
|
|
||||||
if last_ele_signum != ele_signum {
|
|
||||||
// Do 50 iterations of newton's method, should be more than accurate
|
|
||||||
let x = {
|
|
||||||
let mut x1: f64 = last_ele.unwrap().x;
|
|
||||||
let mut x2: f64;
|
|
||||||
let mut fail: bool = false;
|
|
||||||
loop {
|
|
||||||
x2 = x1
|
|
||||||
- (self.function.get_derivative_1(x1)
|
|
||||||
/ self.function.get_derivative_2(x1));
|
|
||||||
if !(self.min_x..self.max_x).contains(&x2) {
|
|
||||||
fail = true;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
if (x2 - x1).abs() < resolution {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
|
|
||||||
x1 = x2;
|
|
||||||
}
|
|
||||||
|
|
||||||
match fail {
|
|
||||||
true => f64::NAN,
|
|
||||||
false => x1,
|
|
||||||
}
|
|
||||||
};
|
|
||||||
|
|
||||||
if !x.is_nan() {
|
|
||||||
extrama_list.push(Value::new(x, self.function.get(x)));
|
|
||||||
}
|
|
||||||
}
|
|
||||||
last_ele = Some(*ele);
|
|
||||||
}
|
|
||||||
self.output.extrema = Some(extrama_list);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
pub fn display(
|
pub fn display(
|
||||||
|
|||||||
67
src/misc.rs
67
src/misc.rs
@ -1,3 +1,7 @@
|
|||||||
|
use std::ops::Range;
|
||||||
|
|
||||||
|
use eframe::egui::plot::Value;
|
||||||
|
|
||||||
cfg_if::cfg_if! {
|
cfg_if::cfg_if! {
|
||||||
if #[cfg(target_arch = "wasm32")] {
|
if #[cfg(target_arch = "wasm32")] {
|
||||||
use wasm_bindgen::prelude::*;
|
use wasm_bindgen::prelude::*;
|
||||||
@ -90,3 +94,66 @@ pub fn digits_precision(x: f64, digits: usize) -> f64 {
|
|||||||
let large_number: f64 = 10.0_f64.powf(digits as f64);
|
let large_number: f64 = 10.0_f64.powf(digits as f64);
|
||||||
(x * large_number).round() / large_number
|
(x * large_number).round() / large_number
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
|
/// Implements newton's method of finding roots.
|
||||||
|
/// `threshold` is the target accuracy threshold
|
||||||
|
/// `range` is the range of valid x values (used to stop calculation when the point won't display anyways)
|
||||||
|
/// `data` is the data to iterate over (a Vector of egui's `Value` struct)
|
||||||
|
/// `f` is f(x)
|
||||||
|
/// `f_` is f'(x)
|
||||||
|
/// The function returns a list of `x` values where roots occur
|
||||||
|
pub fn newtons_method(
|
||||||
|
threshold: f64, range: Range<f64>, data: Vec<Value>, f: &dyn Fn(f64) -> f64,
|
||||||
|
f_1: &dyn Fn(f64) -> f64,
|
||||||
|
) -> Vec<f64> {
|
||||||
|
let mut output_list: Vec<f64> = Vec::new();
|
||||||
|
let mut last_ele: Option<Value> = None;
|
||||||
|
for ele in data.iter() {
|
||||||
|
if last_ele.is_none() {
|
||||||
|
last_ele = Some(*ele);
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
let last_ele_signum = last_ele.unwrap().y.signum();
|
||||||
|
let ele_signum = ele.y.signum();
|
||||||
|
|
||||||
|
// If either are NaN, just continue iterating
|
||||||
|
if last_ele_signum.is_nan() | ele_signum.is_nan() {
|
||||||
|
continue;
|
||||||
|
}
|
||||||
|
|
||||||
|
if last_ele_signum != ele_signum {
|
||||||
|
// Do 50 iterations of newton's method, should be more than accurate
|
||||||
|
let x = {
|
||||||
|
let mut x1: f64 = last_ele.unwrap().x;
|
||||||
|
let mut x2: f64;
|
||||||
|
let mut fail: bool = false;
|
||||||
|
loop {
|
||||||
|
x2 = x1 - (f(x1) / f_1(x1));
|
||||||
|
if !range.contains(&x2) {
|
||||||
|
fail = true;
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
if (x2 - x1).abs() < threshold {
|
||||||
|
break;
|
||||||
|
}
|
||||||
|
|
||||||
|
x1 = x2;
|
||||||
|
}
|
||||||
|
|
||||||
|
match fail {
|
||||||
|
true => f64::NAN,
|
||||||
|
false => x1,
|
||||||
|
}
|
||||||
|
};
|
||||||
|
|
||||||
|
if !x.is_nan() {
|
||||||
|
output_list.push(x);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
last_ele = Some(*ele);
|
||||||
|
}
|
||||||
|
output_list
|
||||||
|
}
|
||||||
|
|||||||
Loading…
x
Reference in New Issue
Block a user