multicalc

Rust 科学计算，用于高精度单变量和多变量微积分

主要特性

• 使用纯 Rust 编写，安全可靠
• 无标准库，零堆分配，无 panic
• 基于特质的泛型实现，支持浮点数和复数
• 完全文档化，包括代码示例
• 全面测试套件，涵盖所有可能的错误条件
• 支持线性、多项式、三角、指数和任何复杂方程，包括任意数量的变量！
  • 任意阶数的高阶数值微分
    • 有限差分法，用于总微分和偏微分
  • 任意阶数的数值积分
    • 迭代方法：Booles、Simpsons、梯形法
    • 高斯求积：高斯-勒让德、高斯-埃尔米特、高斯-拉盖尔
  • 雅可比矩阵和海森矩阵
  • 矢量场微积分：线积分和流量积分，旋度和散度
  • 将任何给定方程逼近线性或二次模式

目录

• 1. 单变量方程的导数
• 2. 多变量方程的导数
• 3. 单变量方程的积分
• 4. 多变量方程的积分
• 5. 雅可比矩阵
• 6. 海森矩阵
• 7. 线性逼近
• 8. 二次逼近
• 9. 线积分和流量积分
• 10. 旋度和散度
• 11. 错误处理
• 12. 试验性

1. 单变量方程的导数

// assume we want to differentiate f(x) = x^3
let func = | arg: f64 | -> f64 
{ 
    return arg*arg*arg;
};

let point = 2.0; //the point at which we want to differentiate

use multicalc::numerical_derivative::derivator::*;
use multicalc::numerical_derivative::finite_difference::*;
 
let derivator = SingleVariableSolver::default();
//alternatively, you can also create the derivator with custom parameters using SingleVariableSolver::from_parameters()

let val = derivator.get(1, &my_func, point).unwrap(); //single derivative
assert!(f64::abs(val - 12.0) < 1e-7);
let val = derivator.get(2, &my_func, point).unwrap(); //double derivative
assert!(f64::abs(val - 12.0) < 1e-5);
let val = derivator.get(3, &my_func, point).unwrap(); //triple derivative
assert!(f64::abs(val - 6.0) < 1e-3);

//for single and double derivatives, you can also just use the convenience wrappers
let val = derivator.get_single(&my_func, point).unwrap(); //single derivative
assert!(f64::abs(val - 12.0) < 1e-7);
let val = derivator.get_double(&my_func, point).unwrap(); //double derivative
assert!(f64::abs(val - 12.0) < 1e-5);

2. 多变量方程的导数

//function is f(x,y,z) = y*sin(x) + x*cos(y) + x*y*e^z
let func = | args: &[f64; 3] | -> f64 
{ 
    return args[1]*args[0].sin() + args[0]*args[1].cos() + args[0]*args[1]*args[2].exp();
};

let point = [1.0, 2.0, 3.0]; //the point of interest

use multicalc::numerical_derivative::derivator::*;
use multicalc::numerical_derivative::finite_difference::*;

let derivator = MultiVariableSolver::default();
//alternatively, you can also create the derivator with custom parameters using MultiVariableSolver::from_parameters()

let idx: [usize; 2] = [0, 1]; //mixed partial double derivate d(df/dx)/dy
let val = derivator.get(2, &my_func, &idx, &point).unwrap();
let expected_value = f64::cos(1.0) - f64::sin(2.0) + f64::exp(3.0);
assert!(f64::abs(val - expected_value) < 0.001);
 
let idx: [usize; 2] = [1, 1]; //partial double derivative d(df/dy)/dy 
let val = derivator.get(2, &my_func, &idx, &point).unwrap();
let expected_value = -1.0*f64::cos(2.0);
assert!(f64::abs(val - expected_value) < 0.0001);

3. 单变量方程的积分

//integrate 2*x . the function would be:
let my_func = | arg: f64 | -> f64 
{ 
    return 2.0*arg;
};

use multicalc::numerical_integration::integrator::*;
use multicalc::numerical_integration::iterative_integration;

let integrator = iterative_integration::SingleVariableSolver::default();  
//alternatively, create a custom integrator using SingleVariableSolver::from_parameters()

//instead of using iterative algorithms, you can also use gaussian quadratures instead:
//let integrator = gaussian_integration::SingleVariableSolver::default();  

let integration_limit = [[0.0, 2.0]; 1]; //desired integration limit 
let val = integrator.get(1, &my_func, &integration_limit).unwrap(); //single integration
assert!(f64::abs(val - 4.0) < 1e-6);

let integration_limit = [[0.0, 2.0], [-1.0, 1.0]]; //desired integration limits
let val = integrator.get(2, &my_func, &integration_limit).unwrap(); //double integration
assert!(f64::abs(val - 8.0) < 1e-6);
///
let integration_limit = [[0.0, 2.0], [0.0, 2.0], [0.0, 2.0]]; //desired integration limits
let val = integrator.get(3, &my_func, &integration_limit).unwrap(); //triple integration
assert!(f64::abs(val - 16.0) < 1e-6);


//you can also just use convenience wrappers for single and double integration
let integration_limit = [0.0, 2.0];
let val = integrator.get_single(&my_func, &integration_limit).unwrap(); //single integration
assert!(f64::abs(val - 4.0) < 1e-6);

let integration_limit = [[0.0, 2.0], [-1.0, 1.0]];
let val = integrator.get_double(&my_func, &integration_limit).unwrap(); //double integration
assert!(f64::abs(val - 8.0) < 1e-6);

4. 多变量方程的积分

//equation is 2.0*x + y*z
let func = | args: &[f64; 3] | -> f64 
{ 
    return 2.0*args[0] + args[1]*args[2];
};

//for multivariable integration to work, we must know the final value of all variables
//for variables that are constant, they will have their constant values
//for variables that are being integrated, their values will be the final upper limit of integration
//in this example, we are keeping 'z' constant at 3.0, and integrating x and y to upper limits 1.0 and 2.0 respectively
let point = [1.0, 2.0, 3.0];

let iterator = iterative_integration::MultiVariableSolver::default();
//you can also create a custom integrator using MultiVariableSolver::from_parameters()

//instead of using iterative algorithms, you can also use gaussian quadratures instead:
//let integrator = gaussian_integration::MultiVariableSolver::default();  

//integration for x, known to be x*x + x*y*z, expect a value of ~7.00
let integration_limit = [0.0, 1.0];
let val = iterator.get_single_partial(&func, 0, &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(1, [0], &func, &[integration_limit; 1], &point);
assert!(f64::abs(val - 7.0) < 1e-7);


//integration for y, known to be 2.0*x*y + y*y*z/2.0, expect a value of ~10.00 
let integration_limit = [0.0, 2.0];
let val = iterator.get_single_partial(&func, 1, &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(1, [1], &func, &[integration_limit; 1], &point);
assert!(f64::abs(val - 10.0) < 1e-7);


//double integration for first x then y, expect a value of ~8.0
let integration_limit = [[0.0, 1.0], [0.0, 2.0]];
let val = integrator.get_double_partial(&func, [0, 1], &integration_limit, &point).unwrap();
//alternative syntax: iterator.get(2, [0, 1], &func, &integration_limit, &point);
assert!(f64::abs(val - 8.0) < 1e-7);

//triple integration for x, expect a value of ~7.0
let integration_limit = [[0.0, 1.0], [0.0, 1.0], [0.0, 1.0]];
let val = integrator.get(3, [0, 0, 0], &func, &integration_limit, &point).unwrap();
assert!(f64::abs(val - 7.0) < 1e-7);

5. 雅可比矩阵

//function is x*y*z
let func1 = | args: &[f64; 3] | -> f64 
{ 
    return args[0]*args[1]*args[2];
};

//function is x^2 + y^2
let func2 = | args: &[f64; 3] | -> f64 
{ 
    return args[0].powf(2.0) + args[1].powf(2.0);
};

let function_matrix: [&dyn Fn(&[f64; 3]) -> f64; 2] = [&func1, &func2];

let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let jacobian = Jacobian::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = jacobian.get(&function_matrix, &points).unwrap();

let expected_result = [[6.0, 3.0, 2.0], [2.0, 4.0, 0.0]];
for i in 0..function_vector.len()
{
    for j in 0..points.len()
    {
        //numerical error less than 1e-6
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 1e-6);
    }
}

6. 海森矩阵

//function is y*sin(x) + 2*x*e^y
let func = | args: &[f64; 2] | -> f64 
{ 
    return args[1]*args[0].sin() + 2.0*args[0]*args[1].exp();
};

let points = [1.0, 2.0]; //the point around which we want the hessian matrix

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let hessian = Hessian::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = hessian.get(&func, &points).unwrap();

assert!(result.len() == points.len()); //number of rows
assert!(result[0].len() == points.len()); //number of columns

let expected_result = [[-2.0*f64::sin(1.0), f64::cos(1.0) + 2.0*f64::exp(2.0)], 
                                        [f64::cos(1.0) + 2.0*f64::exp(2.0), 2.0*f64::exp(2.0)]];

for i in 0..points.len()
{
    for j in 0..points.len()
    {
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 1e-5);
    }
}

7. 线性逼近

//function is x + y^2 + z^3, which we want to linearize
let function_to_approximate = | args: &[f64; 3] | -> f64 
{ 
    return args[0] + args[1].powf(2.0) + args[2].powf(3.0);
};

let point = [1.0, 2.0, 3.0]; //the point we want to linearize around

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let approximator = LinearApproximator::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = approximator.get(&function_to_approximate, &point).unwrap();
assert!(f64::abs(function_to_approximate(&point) - result.get_prediction_value(&point)) < 1e-9);

8. 二次逼近

//function is e^(x/2) + sin(y) + 2.0*z
let function_to_approximate = | args: &[f64; 3] | -> f64 
{ 
    return f64::exp(args[0]/2.0) + f64::sin(args[1]) + 2.0*args[2];
};

let point = [0.0, 3.14/2.0, 10.0]; //the point we want to approximate around

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;

let approximator = QuadraticApproximator::<MultiVariableSolver>::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

let result = approximator.get(&function_to_approximate, &point).unwrap();

assert!(f64::abs(function_to_approximate(&point) - result.get_prediction_value(&point)) < 1e-9);

9. 线积分和流量积分

//vector field is (y, -x). On a 2D plane this would like a tornado rotating counter-clockwise
//curve is a unit circle, defined by (Cos(t), Sin(t))
//limit t goes from 0->2*pi

let vector_field_matrix: [&dyn Fn(&f64, &f64) -> f64; 2] = [&(|_:&f64, y:&f64|-> f64 { *y }), &(|x:&f64, _:&f64|-> f64 { -x })];

let transformation_matrix: [&dyn Fn(&f64) -> f64; 2] = [&(|t:&f64|->f64 { t.cos() }), &(|t:&f64|->f64 { t.sin() })];

let integration_limit = [0.0, 6.28];

//line integral of a unit circle curve on our vector field from 0 to 2*pi, expect an answer of -2.0*pi
let val = line_integral::get_2d(&vector_field_matrix, &transformation_matrix, &integration_limit).unwrap();
assert!(f64::abs(val + 6.28) < 0.01);

//flux integral of a unit circle curve on our vector field from 0 to 2*pi, expect an answer of 0.0
let val = flux_integral::get_2d(&vector_field_matrix, &transformation_matrix, &integration_limit).unwrap();
assert!(f64::abs(val - 0.0) < 0.01);

10. 旋度和散度

//vector field is (2*x*y, 3*cos(y))
//x-component
let vf_x = | args: &[f64; 2] | -> f64 
{ 
    return 2.0*args[0]*args[1];
};

//y-component
let vf_y = | args: &[f64; 2] | -> f64 
{ 
    return 3.0*args[1].cos()
};

let vector_field_matrix: [&dyn Fn(&[f64; 2]) -> f64; 2] = [&vf_x, &vf_y];

let point = [1.0, 3.14]; //the point of interest

use multicalc::numerical_derivative::finite_difference::MultiVariableSolver;
let derivator = MultiVariableSolver::default();
//instead of using the in-built MultiVariableSolver, you can also supply your own by implementing the base trait

//curl is known to be -2*x, expect and answer of -2.0
let val = curl::get_2d(derivator, &vector_field_matrix, &point).unwrap();
assert!(f64::abs(val + 2.0) < 0.000001); //numerical error less than 1e-6

//divergence is known to be 2*y - 3*sin(y), expect and answer of 6.27
let val = divergence::get_2d(derivator, &vector_field_matrix, &point).unwrap();
assert!(f64::abs(val - 6.27) < 0.01);

11. 错误处理

尽可能提供“安全”版本的函数，填充默认值并直接返回所需的解。然而，这并不总是可能的，因为没有默认参数可以假设，或者对于故意给用户调整参数自由的函数。在这种情况下，返回一个 Result<T, &'static str> 对象，其中可以查看所有可能的 &'static str 在 错误代码。

12. 试验性

启用“堆”功能以访问某些模块中的基于 std::Vec 的方法。目前，这仅支持通过 get_on_heap() 和 get_on_heap_custom() 方法访问的 雅可比 模块。输出是一个动态分配的 Vec<Vec<T>>。这是为了支持大型数据集，否则可能因静态数组而导致堆栈溢出。未来的计划可能包括为逼近模块添加此类支持。

//function is x*y*z
let func1 = | args: &[f64; 3] | -> f64 
{ 
    return args[0]*args[1]*args[2];
};

//function is x^2 + y^2
let func2 = | args: &[f64; 3] | -> f64 
{ 
    return args[0].powf(2.0) + args[1].powf(2.0);
};

let function_matrix: Vec<Box<dyn Fn(&[f64; 3]) -> f64>> = std::vec![Box::new(func1), Box::new(func2)];

let points = [1.0, 2.0, 3.0]; //the point around which we want the jacobian matrix

let jacobian = Jacobian::<MultiVariableSolver>::default();

let result: Vec<Vec<f64>> = jacobian.get_on_heap(&function_matrix, &points).unwrap();

assert!(result.len() == function_matrix.len()); //number of rows
assert!(result[0].len() == points.len()); //number of columns

let expected_result = [[6.0, 3.0, 2.0], [2.0, 4.0, 0.0]];

for i in 0..function_matrix.len()
{
    for j in 0..points.len()
    {
        assert!(f64::abs(result[i][j] - expected_result[i][j]) < 0.000001);
    }
}

基准测试

见 BENCHMARKS.md

贡献指南

见 CONTRIBUTIONS.md

许可证

multicalc 在 MIT 许可证下授权。

致谢

multicalc 使用 num-complex 为所有浮点类型和复数提供通用功能

联系方式

[email protected]

路线图

• 添加高斯-埃尔米特
• 添加高斯-拉盖尔
• [] 添加用户友好的宏以简化使用
• [] 将无限积分限制添加到迭代积分方法中
• [] 将有限积分限制添加到高斯-埃尔米特
• [] 将有限积分限制添加到高斯-拉盖尔
• [] 将复数支持添加到积分模块中
• [] 添加常微分方程求解器模块
• [] 重构向量微积分模块