C#Unity3D微积分

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

只是想问一下这些函数在数值上计算导数和偏导数是否正确,以及如何像在薛定inger方程,热方程等公式中那样进行二阶偏导和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);
  }
}

如果您经常使用衍生工具,则可以尝试将其声明为扩展方法：

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在Sin函数的[0 .. 2 * PI)范围内时的最大误差

// 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

例如,x = 0时的d(abs(x))/ dx

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处计算导数时,x处不需要f.例如

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

请注意,存在x = 0处的导数df / dx(它是正无穷大).这就是为什么在计算导数时应避免计算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)

相关推荐