1. 数学解析与C++的完美结合
作为一名长期从事科学计算开发的工程师，我深知数学解析在工程实践中的重要性。最近因为腿伤休息了一段时间，现在终于可以重新投入工作，迫不及待想和大家分享这个完整的C++数学解析实现方案。
数学解析（Mathematical Analysis）是研究连续量变化规律的数学分支，它构成了现代科学计算的基石。从物理引擎到金融建模，从信号处理到机器学习，几乎每个工程领域都离不开这些数学工具。而C++以其卓越的性能和灵活性，成为实现这些数学算法的首选语言。
2. C++数学基础库搭建
2.1 核心头文件与常量定义
任何数学计算项目都需要从基础搭建开始。以下是我们的核心头文件包含和常量定义：
cpp复制#include <iostream>
#include <vector>
#include <cmath>
#include <iomanip>
#include <functional>
#include <complex>
#include <algorithm>
#include <stdexcept>
using namespace std;
// 数学常量定义
const double EPS = 1e-12; // 浮点比较精度阈值
const double PI = acos(-1.0); // 圆周率
const double E = exp(1.0); // 自然对数底
// 常用类型别名
using Func = function<double(double)>;
using Vec = vector<double>;
using Mat = vector<Vec>;
using Complex = complex<double>;
提示：EPS值的设置需要根据具体应用场景调整。对于大多数工程计算，1e-12提供了良好的精度平衡。但在极端情况下，可能需要减小到1e-15或增大到1e-8。
2.2 高精度计算注意事项
在实现数学算法时，浮点精度问题是我们需要特别注意的：
- 避免大数吃小数：当两个数量级相差很大的数相加时，较小的数可能会被忽略
- 警惕累积误差：在迭代计算中，误差会不断累积
- 合理选择步长：数值微分和积分中的步长h需要仔细选择
cpp复制// 安全的浮点数比较函数
bool almost_equal(double a, double b, double eps = EPS) {
return fabs(a - b) < eps;
}
3. 极限理论与数值计算
3.1 极限的数值逼近
数学分析中的极限概念可以通过数值方法进行近似计算。我们采用双侧逼近法：
cpp复制double limit(Func f, double a, double h = 1e-6) {
double left = f(a - h);
double right = f(a + h);
return (left + right) * 0.5;
}
这个实现简单但有效，对于大多数连续函数都能给出合理的结果。测试几个经典极限：
cpp复制// 测试重要极限
void test_limits() {
auto lim1 = limit([](double x) { return sin(x)/x; }, 0);
cout << "lim(sin(x)/x) as x->0: " << lim1 << endl;
auto lim2 = important_limit_e();
cout << "lim(1+1/x)^x as x->∞: " << lim2 << endl;
}
3.2 重要极限的特殊实现
对于某些已知形式的极限，我们可以采用更高效的特化实现：
cpp复制double important_limit_e(int n = 1000000) {
return pow(1.0 + 1.0 / n, n);
}
这个实现直接计算了自然对数底e的近似值，通过增大n可以提高精度。
4. 单变量微积分实现
4.1 数值微分技术
导数的中心差分法是工程中最常用的数值微分方法：
cpp复制double derivative(Func f, double x, double h = 1e-5) {
return (f(x + h) - f(x - h)) / (2 * h);
}
步长h的选择至关重要：
- h太大：截断误差大
- h太小：舍入误差大
- 理想h：约为√(EPS)*x
4.2 高阶导数计算
通过递归我们可以实现任意阶导数的计算：
cpp复制double nth_derivative(Func f, double x, int n, double h = 1e-4) {
if (n == 0) return f(x);
if (n == 1) return derivative(f, x, h);
return (nth_derivative(f, x+h, n-1, h) -
nth_derivative(f, x-h, n-1, h)) / (2*h);
}
注意：高阶导数的数值计算精度会迅速下降，通常不建议超过4阶。
5. 多变量微积分实现
5.1 偏导数与梯度
对于多变量函数，我们首先实现偏导数计算：
cpp复制double partial_derivative(function<double(Vec)> f, Vec x, int dim, double h = 1e-5) {
Vec x1 = x, x2 = x;
x1[dim] += h;
x2[dim] -= h;
return (f(x1) - f(x2)) / (2 * h);
}
基于偏导数，梯度计算就水到渠成了：
cpp复制Vec gradient(function<double(Vec)> f, Vec x, double h = 1e-5) {
Vec grad(x.size());
for (int i = 0; i < x.size(); i++) {
grad[i] = partial_derivative(f, x, i, h);
}
return grad;
}
5.2 雅可比矩阵与海森矩阵
雅可比矩阵是多变量函数的导数矩阵，实现如下：
cpp复制Mat jacobian(function<Vec(Vec)> f, Vec x, double h = 1e-5) {
int m = f(x).size();
int n = x.size();
Mat J(m, Vec(n));
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
Vec x1 = x, x2 = x;
x1[j] += h;
x2[j] -= h;
J[i][j] = (f(x1)[i] - f(x2)[i]) / (2*h);
}
}
return J;
}
海森矩阵是二阶偏导数的对称矩阵，可用于优化算法：
cpp复制Mat hessian(function<double(Vec)> f, Vec x, double h = 1e-5) {
int n = x.size();
Mat H(n, Vec(n));
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
Vec x1 = x, x2 = x, x3 = x, x4 = x;
x1[i] += h; x1[j] += h;
x2[i] += h; x2[j] -= h;
x3[i] -= h; x3[j] += h;
x4[i] -= h; x4[j] -= h;
H[i][j] = (f(x1) - f(x2) - f(x3) + f(x4)) / (4*h*h);
}
}
return H;
}
6. 数值积分方法比较
6.1 矩形法与梯形法
最简单的数值积分方法是矩形法：
cpp复制double integrate_rect(Func f, double a, double b, int n = 10000) {
double h = (b - a) / n;
double sum = 0;
for (int i = 0; i < n; i++) sum += f(a + i*h);
return sum * h;
}
梯形法是对矩形法的改进：
cpp复制double integrate_trap(Func f, double a, double b, int n = 10000) {
double h = (b - a) / n;
double sum = 0.5 * (f(a) + f(b));
for (int i = 1; i < n; i++) sum += f(a + i*h);
return sum * h;
}
6.2 辛普森积分法
辛普森法利用抛物线近似，精度更高：
cpp复制double integrate_simpson(Func f, double a, double b, int n = 10000) {
if (n % 2 == 1) n++;
double h = (b - a) / n;
double sum = f(a) + f(b);
for (int i = 1; i < n; i++) {
sum += f(a + i*h) * (i%2 ? 4 : 2);
}
return sum * h / 3.0;
}
6.3 自适应积分策略
对于变化剧烈的函数，自适应积分更高效：
cpp复制double integrate_adaptive(Func f, double a, double b, double tol = 1e-6) {
double h = (b - a) / 2;
double fa = f(a), fb = f(b), fm = f(a + h);
double S = (fa + 4*fm + fb) * h / 3;
double fl = f(a + h/2), fr = f(a + 3*h/2);
double S2 = (fa + 4*fl + 2*fm + 4*fr + fb) * h / 6;
if (fabs(S - S2) < tol) return S2;
return integrate_adaptive(f, a, a+h, tol/2) +
integrate_adaptive(f, a+h, b, tol/2);
}
7. 泰勒级数与函数逼近
泰勒级数是函数在某点附近的幂级数展开：
cpp复制double taylor_series(Func f, double x, double a, int order = 5) {
double res = 0;
double fact = 1;
for (int n = 0; n <= order; n++) {
if (n > 0) fact *= n;
double d = nth_derivative(f, a, n);
res += d * pow(x - a, n) / fact;
}
return res;
}
实际应用中，我们可以用泰勒级数实现各种数学函数：
cpp复制double my_sin(double x) {
// 将x转换到[-π, π]区间
x = fmod(x, 2*PI);
if (x > PI) x -= 2*PI;
else if (x < -PI) x += 2*PI;
double sum = 0;
double term = x;
int n = 1;
while (fabs(term) > EPS) {
sum += term;
term *= -x*x / ((2*n) * (2*n+1));
n++;
}
return sum;
}
8. 常微分方程数值解法
8.1 欧拉方法
最简单的ODE解法是欧拉方法：
cpp复制Vec ode_euler(Func f, double y0, double t0, double t1, int steps) {
Vec res;
double h = (t1 - t0) / steps;
double y = y0, t = t0;
res.push_back(y);
for (int i = 0; i < steps; i++) {
y += h * f(t);
t += h;
res.push_back(y);
}
return res;
}
8.2 龙格-库塔法(RK4)
RK4方法精度更高，是工程实践中的标准选择：
cpp复制double rk4_step(function<double(double,double)> f, double t, double y, double h) {
double k1 = f(t, y);
double k2 = f(t + h/2, y + h*k1/2);
double k3 = f(t + h/2, y + h*k2/2);
double k4 = f(t + h, y + h*k3);
return y + h/6*(k1 + 2*k2 + 2*k3 + k4);
}
Vec ode_rk4(function<double(double,double)> f, double y0, double t0, double t1, int steps) {
Vec res;
double h = (t1 - t0) / steps;
double y = y0, t = t0;
res.push_back(y);
for (int i = 0; i < steps; i++) {
y = rk4_step(f, t, y, h);
t += h;
res.push_back(y);
}
return res;
}
9. 傅里叶分析与信号处理
9.1 离散傅里叶变换(DFT)
DFT是信号处理的基础工具：
cpp复制vector<Complex> dft(const vector<Complex>& a) {
int n = a.size();
vector<Complex> res(n);
for (int k = 0; k < n; k++) {
Complex sum = 0;
for (int t = 0; t < n; t++) {
double angle = -2 * PI * k * t / n;
sum += a[t] * Complex(cos(angle), sin(angle));
}
res[k] = sum;
}
return res;
}
9.2 快速傅里叶变换(FFT)
FFT是DFT的高效实现：
cpp复制vector<Complex> fft(const vector<Complex>& a) {
int n = a.size();
if (n == 1) return a;
vector<Complex> even(n/2), odd(n/2);
for (int i = 0; i < n/2; i++) {
even[i] = a[2*i];
odd[i] = a[2*i+1];
}
auto q = fft(even);
auto r = fft(odd);
vector<Complex> y(n);
for (int k = 0; k < n/2; k++) {
double ang = -2 * PI * k / n;
Complex wk = Complex(cos(ang), sin(ang));
y[k] = q[k] + wk * r[k];
y[k + n/2] = q[k] - wk * r[k];
}
return y;
}
10. 最优化算法实现
10.1 梯度下降法
最基本的优化算法是梯度下降：
cpp复制Vec gradient_descent(function<double(Vec)> f, Vec x0, double lr = 0.01, int iter = 1000) {
Vec x = x0;
for (int i = 0; i < iter; i++) {
Vec g = gradient(f, x);
for (int j = 0; j < x.size(); j++) {
x[j] -= lr * g[j];
}
}
return x;
}
10.2 牛顿法
牛顿法利用二阶导数信息，收敛更快：
cpp复制double newton_method(Func f, double x0, int max_iter = 100) {
double x = x0;
for (int i = 0; i < max_iter; i++) {
double d = derivative(f, x);
if (fabs(d) < EPS) break;
double d2 = nth_derivative(f, x, 2);
x -= d / d2;
}
return x;
}
11. 数学解析的工程应用
11.1 物理引擎开发
刚体动力学本质上就是微分方程的求解问题。以简单的弹簧系统为例：
cpp复制void simulate_spring() {
double k = 1.0; // 弹性系数
double m = 1.0; // 质量
double b = 0.1; // 阻尼系数
auto spring_eq = [k, m, b](double t, pair<double,double> state) {
double x = state.first;
double v = state.second;
return make_pair(v, (-k*x - b*v)/m);
};
// 使用RK4方法求解
pair<double,double> state(1.0, 0.0); // 初始位移和速度
double t = 0, dt = 0.01;
for (int i = 0; i < 1000; i++) {
auto k1 = spring_eq(t, state);
auto k2 = spring_eq(t + dt/2, {
state.first + dt/2 * k1.first,
state.second + dt/2 * k1.second
});
auto k3 = spring_eq(t + dt/2, {
state.first + dt/2 * k2.first,
state.second + dt/2 * k2.second
});
auto k4 = spring_eq(t + dt, {
state.first + dt * k3.first,
state.second + dt * k3.second
});
state.first += dt/6 * (k1.first + 2*k2.first + 2*k3.first + k4.first);
state.second += dt/6 * (k1.second + 2*k2.second + 2*k3.second + k4.second);
t += dt;
cout << t << " " << state.first << endl;
}
}
11.2 机器学习中的反向传播
神经网络训练的核心是梯度计算，这正是多变量微积分的应用：
cpp复制class NeuralNetwork {
vector<Mat> weights;
vector<Vec> biases;
Vec forward(const Vec& input) {
Vec activation = input;
for (size_t i = 0; i < weights.size(); i++) {
Vec z(weights[i].size());
for (size_t j = 0; j < weights[i].size(); j++) {
z[j] = inner_product(weights[i][j].begin(), weights[i][j].end(),
activation.begin(), biases[i][j]);
z[j] = 1.0 / (1.0 + exp(-z[j])); // sigmoid
}
activation = z;
}
return activation;
}
pair<vector<Mat>, vector<Vec>> backward(const Vec& input, const Vec& target) {
// 前向传播
vector<Vec> activations = {input};
vector<Vec> zs;
Vec activation = input;
for (size_t i = 0; i < weights.size(); i++) {
Vec z(weights[i].size());
for (size_t j = 0; j < weights[i].size(); j++) {
z[j] = inner_product(weights[i][j].begin(), weights[i][j].end(),
activation.begin(), biases[i][j]);
}
zs.push_back(z);
activation.resize(z.size());
for (size_t j = 0; j < z.size(); j++) {
activation[j] = 1.0 / (1.0 + exp(-z[j]));
}
activations.push_back(activation);
}
// 反向传播
vector<Mat> nabla_w(weights.size());
vector<Vec> nabla_b(biases.size());
Vec delta = activations.back();
for (size_t j = 0; j < delta.size(); j++) {
delta[j] = (delta[j] - target[j]) * delta[j] * (1 - delta[j]);
}
for (int i = weights.size()-1; i >= 0; i--) {
nabla_b[i] = delta;
nabla_w[i].resize(weights[i].size());
for (size_t j = 0; j < weights[i].size(); j++) {
nabla_w[i][j].resize(weights[i][j].size());
for (size_t k = 0; k < weights[i][j].size(); k++) {
nabla_w[i][j][k] = activations[i][k] * delta[j];
}
}
if (i > 0) {
Vec new_delta(activations[i].size());
for (size_t k = 0; k < new_delta.size(); k++) {
double sum = 0;
for (size_t j = 0; j < delta.size(); j++) {
sum += weights[i][j][k] * delta[j];
}
new_delta[k] = sum * activations[i][k] * (1 - activations[i][k]);
}
delta = new_delta;
}
}
return {nabla_w, nabla_b};
}
};
12. 性能优化与工程实践
12.1 模板元编程优化
对于性能关键的数学计算，我们可以使用模板元编程技术：
cpp复制template <typename T, int N>
class Vector {
T data[N];
public:
Vector() { fill(data, data+N, T()); }
template <typename... Args>
Vector(Args... args) : data{static_cast<T>(args)...} {}
T& operator[](int i) { return data[i]; }
const T& operator[](int i) const { return data[i]; }
Vector operator+(const Vector& other) const {
Vector result;
for (int i = 0; i < N; i++) {
result[i] = data[i] + other[i];
}
return result;
}
// 其他运算符重载...
};
// 使用示例
Vector<double, 3> v1(1.0, 2.0, 3.0);
Vector<double, 3> v2(4.0, 5.0, 6.0);
auto v3 = v1 + v2;
12.2 SIMD并行化
现代CPU支持SIMD指令，可以大幅提升数学计算性能：
cpp复制#include <immintrin.h>
void vector_add(float* a, float* b, float* c, int n) {
for (int i = 0; i < n; i += 8) {
__m256 va = _mm256_load_ps(a + i);
__m256 vb = _mm256_load_ps(b + i);
__m256 vc = _mm256_add_ps(va, vb);
_mm256_store_ps(c + i, vc);
}
}
13. 测试与验证策略
13.1 单元测试框架
为数学库编写全面的测试用例至关重要：
cpp复制#define ASSERT_NEAR(a, b, eps) \
if (fabs((a)-(b)) > (eps)) { \
cerr << "Assertion failed: " << #a << " (" << (a) << ") != " << #b << " (" << (b) << ")" << endl; \
return false; \
}
bool test_derivative() {
auto f = [](double x) { return x*x; };
double d = derivative(f, 2.0);
ASSERT_NEAR(d, 4.0, 1e-5);
return true;
}
bool test_integral() {
auto f = [](double x) { return x*x; };
double i = integrate_simpson(f, 0, 2);
ASSERT_NEAR(i, 8.0/3, 1e-5);
return true;
}
void run_tests() {
cout << "Running tests..." << endl;
bool pass = true;
pass &= test_derivative();
pass &= test_integral();
// 更多测试...
if (pass) {
cout << "All tests passed!" << endl;
} else {
cout << "Some tests failed!" << endl;
}
}
13.2 数值稳定性分析
数学算法的数值稳定性需要特别关注：
cpp复制void analyze_stability() {
auto f = [](double x) { return exp(x) - 1; };
cout << "Derivative stability analysis:" << endl;
for (double h = 1e-1; h > 1e-20; h *= 0.1) {
double d = derivative(f, 0, h);
cout << "h = " << h << ", df/dx = " << d << ", error = " << fabs(d - 1.0) << endl;
}
cout << "\nIntegral stability analysis:" << endl;
for (int n = 10; n <= 1e6; n *= 10) {
double i = integrate_simpson(f, 0, 1, n);
cout << "n = " << n << ", integral = " << i << ", error = " << fabs(i - (E-1)) << endl;
}
}
14. 跨平台兼容性考虑
14.1 编译器兼容性处理
不同编译器对数学标准的支持可能有差异：
cpp复制// 确保M_PI等常量定义存在
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif
// 处理不同平台的复杂数支持
#ifdef _MSC_VER
using Complex = _Dcomplex;
#else
using Complex = complex<double>;
#endif
14.2 精度控制与可移植性
cpp复制// 可移植的精度控制函数
template <typename T>
typename enable_if<is_floating_point<T>::value, T>::type
epsilon() {
return numeric_limits<T>::epsilon();
}
// 跨平台的随机数生成
double random_double(double min = 0.0, double max = 1.0) {
static random_device rd;
static mt19937 gen(rd());
uniform_real_distribution<> dis(min, max);
return dis(gen);
}
15. 数学库的扩展与定制
15.1 添加新数学函数
扩展库功能时保持一致的接口风格：
cpp复制namespace special_functions {
// 伽马函数近似
double gamma(double x) {
// Lanczos近似实现
const double g = 7;
static const double p[] = {
0.99999999999980993, 676.5203681218851, -1259.1392167224028,
771.32342877765313, -176.61502916214059, 12.507343278686905,
-0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
};
if (x < 0.5) {
return PI / (sin(PI * x) * gamma(1 - x));
}
x -= 1;
double a = p[0];
double t = x + g + 0.5;
for (int i = 1; i < sizeof(p)/sizeof(p[0]); ++i) {
a += p[i] / (x + i);
}
return sqrt(2*PI) * pow(t, x+0.5) * exp(-t) * a;
}
}
15.2 支持自动微分
对于需要高阶微分的情况，可以实现自动微分：
cpp复制template <typename T>
struct Dual {
T value;
T derivative;
Dual(T v, T d = T()) : value(v), derivative(d) {}
Dual operator+(const Dual& other) const {
return Dual(value + other.value, derivative + other.derivative);
}
Dual operator*(const Dual& other) const {
return Dual(value * other.value,
derivative * other.value + value * other.derivative);
}
// 其他运算符重载...
};
template <typename Func>
auto autodiff(Func f, double x) {
Dual<double> xd(x, 1.0);
auto yd = f(xd);
return yd.derivative;
}
16. 实际工程应用案例
16.1 期权定价模型
金融工程中的Black-Scholes模型实现：
cpp复制double normal_cdf(double x) {
return 0.5 * erfc(-x * M_SQRT1_2);
}
double black_scholes(double S, double K, double T, double r, double sigma, bool is_call) {
double d1 = (log(S/K) + (r + 0.5*sigma*sigma)*T) / (sigma*sqrt(T));
double d2 = d1 - sigma*sqrt(T);
if (is_call) {
return S * normal_cdf(d1) - K * exp(-r*T) * normal_cdf(d2);
} else {
return K * exp(-r*T) * normal_cdf(-d2) - S * normal_cdf(-d1);
}
}
16.2 机器人运动规划
使用微分方程求解机器人轨迹：
cpp复制struct RobotState {
double x, y, theta; // 位置和朝向
double v, omega; // 线速度和角速度
};
RobotState simulate_robot(RobotState state, double dt,
function<pair<double,double>(double,RobotState)> control) {
auto [a, alpha] = control(dt, state);
// 运动学模型
state.v += a * dt;
state.omega += alpha * dt;
state.theta += state.omega * dt;
state.x += state.v * cos(state.theta) * dt;
state.y += state.v * sin(state.theta) * dt;
return state;
}
void plan_trajectory() {
RobotState state{0, 0, 0, 0, 0};
double dt = 0.01;
auto controller = [](double t, RobotState s) {
double a = (t < 1.0) ? 0.5 : -0.5;
double alpha = sin(t);
return make_pair(a, alpha);
};
for (double t = 0; t < 5.0; t += dt) {
state = simulate_robot(state, dt, controller);
cout << t << " " << state.x << " " << state.y << endl;
}
}
17. 性能基准测试
比较不同算法的性能表现：
cpp复制#include <chrono>
void benchmark() {
auto f = [](double x) { return sin(x)*exp(-x*x); };
auto start = chrono::high_resolution_clock::now();
double r1 = integrate_rect(f, -10, 10, 1e6);
auto end = chrono::high_resolution_clock::now();
cout << "Rectangle: " << chrono::duration<double>(end-start).count() << "s, result=" << r1 << endl;
start = chrono::high_resolution_clock::now();
double r2 = integrate_trap(f, -10, 10, 1e6);
end = chrono::high_resolution_clock::now();
cout << "Trapezoid: " << chrono::duration<double>(end-start).count() << "s, result=" << r2 << endl;
start = chrono::high_resolution_clock::now();
double r3 = integrate_simpson(f, -10, 10, 1e6);
end = chrono::high_resolution_clock::now();
cout << "Simpson: " << chrono::duration<double>(end-start).count() << "s, result=" << r3 << endl;
start = chrono::high_resolution_clock::now();
double r4 = integrate_adaptive(f, -10, 10);
end = chrono::high_resolution_clock::now();
cout << "Adaptive: " << chrono::duration<double>(end-start).count() << "s, result=" << r4 << endl;
}
18. 数学可视化技术
18.1 函数绘图工具
将数学函数可视化有助于理解：
cpp复制void plot_function(Func f, double a, double b, int points = 100) {
ofstream data("function.dat");
double dx = (b - a) / points;
for (int i = 0; i <= points; i++) {
double x = a + i * dx;
data << x << " " << f(x) << "\n";
}
data.close();
// 使用Gnuplot绘图
system("gnuplot -p -e \"plot 'function.dat' with lines title 'Function'\"");
}
18.2 向量场可视化
对于多变量函数，可以绘制梯度场：
cpp复制void plot_gradient_field(function<double(Vec)> f,
double xmin, double xmax,
double ymin, double ymax,
int grid = 20) {
ofstream data("gradient.dat");
double dx = (xmax - xmin) / grid;
double dy = (ymax - ymin) / grid;
for (int i = 0; i <= grid; i++) {
for (int j = 0; j <= grid; j++) {
double x = xmin + i * dx;
double y = ymin + j * dy;
Vec grad = gradient(f, {x, y});
// 归一化向量长度以便可视化
double len = sqrt(grad[0]*grad[0] + grad[1]*grad[1]);
if (len > 1e-6) {
grad[0] /= len;
grad[1] /= len;
}
data << x << " " << y << " "
<< grad[0]*dx*0.3 << " " << grad[1]*dy*0.3 << "\n";
}
}
data.close();
system("gnuplot -p -e \""
"set title 'Gradient Field';"
"plot 'gradient.dat' using 1:2:3:4 with vectors head filled lt 2 title 'Gradient'\"");
}
19. 数学符号计算扩展
虽然C++主要是数值计算，但可以简单实现符号计算：
cpp复制class Symbol {
string repr;
function<double(const map<string,double>&)> eval;
public:
Symbol(string r, function<double(const map<string,double>&)> e)
: repr(r), eval(e) {}
double operator()(const map<string,double>& vars) const {
return eval(vars);
}
string str() const { return repr; }
Symbol derivative(const string& var) const {
// 简单实现符号微分
// 实际工程中应使用更复杂的实现
if (repr == var) {
return Symbol("1", [](auto) { return 1.0; });
} else if (repr.find(var) == string::npos) {
return Symbol("0", [](auto) { return 0.0; });
} else {
// 简单处理常见函数
if (repr.find("sin("+var+")") != string::npos) {
return Symbol("cos("+var+")", [](auto vars) {
return cos(vars.at(var));
});
}
// 其他情况简化处理
return Symbol("d("+repr+")/d"+var, [](auto) { return NAN; });
}
}
};
void symbolic_demo() {
Symbol x("x", [](auto vars) { return vars.at("x"); });
Symbol y("y", [](auto vars) { return vars.at("y"); });
auto f = x * y; // 实际应重载运算符实现
cout << "f(x,y) = " << f.str() <<