c# - C#Unity3D微积分

我想使用Unity进行类似于MATLAB的数学运算，而且还希望在化学，物理学和工程学领域使用科学类型的东西。

只是想问一下这些函数对于数字计算导数和偏导数是否正确，以及如何像在薛定2方程，热方程等公式中那样进行二阶偏导和Laplace运算符？

我仍在学习微分方程，但想将其与C＃中的数值计算联系起来进行计算。

public double Derivative(Func<double, double> function, double x, double h)
{
    return (function(x + h) - function(x)) / h;
}

public double SecondDerivative(Func<double, double> function, double x, double h)
{
    return (function(x + h) - 2 * function(x) + function(x - h)) / (h * h);
}

public double ThirdDerivative(Func<double, double> function, double x, double h)
{
    return (function(x + 3 * h) - 3 * function(x + 2 * h) + 3 * function(x + h) - function(x)) / (h * h * h);
}

public double PartialDerivativeX(Func<double, double, double> function, double x, double y, double h)
{
    return (function(x + h, y) - function(x, y)) / h;
}

public double PartialDerivativeY(Func<double, double, double> function, double x, double y, double h)
{
    return (function(x, y + h) - function(x, y)) / h;
}

最佳答案

您的实现是一个很好的近似值（因为导数是一个极限，而h是一个有限值）。但是，我建议使用一些不同的代码：

public static class MyMath {
  // static: we don't want "this"
  // Func<double, double> return value: derivative is a function, not a value.
  //   If we want a point - double - let's name the method as DerivativeAt
  // No h - we can't provide correct h for all possible x
  public static Func<double, double> Derivative(Func<double, double> function) {
    //DONE: Validate public methods arguments
    if (null == function)
      throw new ArgumentNullException("function");

    return new Func<double, double>((x) => {
      // Let's compute h for given x
      // Easiest, but not the best
      double h = Math.Abs(x) < 1e-10 ? 1e-16 : x / 1.0e6;

      // "Central" derivative is often a better choice then right one ((f(x + h) - f(x))/h)
      return (function(x + h) - function(x - h)) / (2.0 * h);
    });
  }

  // h = 0.0: be nice and let user has no idea what step is reasonable
  public static double DerivativeAt(Func<double, double> function,
                                    double x,
                                    double h = 0.0) {
    //DONE: Validate public methods arguments
    if (null == function)
      throw new ArgumentNullException("function");

    // If user don't want to provide h, let's compute it
    if (0 == h)
      h = Math.Abs(x) < 1e-10 ? 1e-16 : x / 1.0e6; // Easiest, but not the best

    // "Central" derivative is often a better choice then right one ((f(x + h) - f(x))/h)
    return (function(x + h) - function(x - h)) / (2.0 * h);
  }
}

如果您经常使用Derivative，则可以尝试将其声明为扩展方法：

public static Func<double, double> Derivative(this Func<double, double> function) {...}

public static double DerivativeAt(this Func<double, double> function,
                                  double x,
                                  double h = 0.0) { ... }

演示：让我们找出x函数的[0 .. 2 * PI)范围在Sin时的最大错误

// We don't want to repeat pesky "MyMath" in "MyMath.Derivative"
using static MyNamespace.MyMath;

...

// Derivative of Sin (expected to be Cos)
var d_sin = Derivative(x => Math.Sin(x));

double maxError = Enumerable
  .Range(0, 1000)
  .Select(i => 2.0 * Math.PI * i / 1000.0)
  .Select(x => Math.Abs(d_sin(x) - Math.Cos(x))) // d(sin(x)) / dx == cos(x)
  .Max();

Console.WriteLine(maxError);

结果：

1.64271596325705E-10

编辑：“中央”派生。

众所周知，导数是一个极限

df/dx == lim (f(x + h) - f(x)) / h
         h -> 0

但是我们可以问：h趋向于0的趋势。在有复数的情况下，我们有很多方法（例如h可以螺旋下降到0或沿着海峡线）；对于实数，h可以为正（右半导数）或负（左半导数）。通常（标准清晰度），我们要求左半导数等于右半导数才能具有导数：

d+f(x) == d-f(x) == df/dx

但是，有时我们使用宽大的定义（“中心”派生）：

df/dx == (d+f(x) + d-f(x)) / 2

例如，在d(abs(x))/dx处的x = 0

d-abs(x)      = -1
d+abs(x)      =  1
d abs(x) / dx    doesn't exist (standard definition)
d abs(x) / dx =  0 "central", lenient definition.

请注意，您当前的代码实际上计算的是右半导数；如果是Abs(x)，您会得到错误的1。如果不是在微积分中，而是在工程学中（想象有速度的运动汽车），0在上下文中是一个更好的答案。另一个问题是，在x处计算导数时，我们不需要在f处存在x。例如

f(x) = x / abs(x) which can be put as

       -1 when x < 0
f(x) =    doesn't exist when x = 0
       +1 when x > 0

请注意，在df/dx处的导数x = 0存在（它是正无穷大）。这就是为什么在计算导数时应避免计算f(x)的原因。您当前的代码将返回

 (h / h + 0 / 0) / h == (1 + NaN) / h == NaN

double.NaN-导数不存在（那是错误的）； “中央”衍生产品将返回

 (h/h - -h/h) / (2 * h) == 1 / h == some huge number (approximation of +Inf)