直观的 GSL 僵硬 EDO

如何解决直观的 GSL 僵硬 EDO

我需要解决一个行为僵硬的 ODE。所述 ODE 将由

给出

int odefunc (double x,const double y[],double f[],void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    double y_eq=-1.931+1.5*log(x)-x;
    f[0] = k*(exp(2*y_eq-y[0])-exp(y[0]));   
    return GSL_SUCCESS;
}

我不太习惯 C（不过我习惯 C++）而且 the GSL page for Differential Equations 中的文档没有太多注释让我理解它。我已经使用 this post 来理解它（也是因为我只有 1 个 ode，而不是它们的系统），并且在某种程度上它很有用，但对我来说仍然有点困惑。对于我的测试，我使用了我的函数而不是给定的函数，尽管我不理解代码，但我已经使用了它。但现在对我来说理解是关键，因为...

我需要修改他们给我的例子，我需要使用一种算法来解决刚性方程，最好是 gsl_odeiv2_step_msbdf 具有自适应步长或等效的东西（似乎使用了 BDF 方法在学术界很常见）。这需要 jacobian，如果你不习惯 GSL 或 C，那部分真的很复杂。

解决方法

要手动实现算法差异化，包含明确的中间步骤会有所帮助

int odefunc (double x,const double y[],double f[],void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    double y_eq=-1.931+1.5*log(x)-x;
    double v1 = exp(2*y_eq-y[0]);
    double v2 = exp(y[0]);
    double z = k*(v1-v2);   
    f[0] = z;
    return GSL_SUCCESS;
}

然后你可以用它的导数传播来增强每个步骤

int jac (double x,double *dfdy,double dfdx[],void *params)
{
    double k=1e12; //At k=1e7 it works (doesn't show stiffness)
    // Dx_dx=Dy_dy=1,Dx_dy=Dy_dx=0 are used without naming them
    double y_eq = -1.931+1.5*log(x)-x,Dy_eq_dx=1.5/x-1,Dy_eq_dy=0;
    double v1 = exp(2*y_eq-y[0]),Dv1_dx=v1*(2*Dy_eq_dx-0),Dv1_dy=v1*(2*Dy_eq_dy-1);
    double v2 = exp(y[0]),Dv2_dx=v2*(0),Dv2_dy=v2*(1);
    double z = k*(v1-v2),Dz_dx=k*(Dv1_dx-Dv2_dx),Dz_dy=k*(Dv1_dy-Dv2_dy);   
    dfdx[0] = Dz_dx;
    dfdy[0] = Dz_dy;
    return GSL_SUCCESS;
}

对于较大的代码，请使用自动执行这些步骤的代码重写工具。 auto-diff.org 应该有一个合适的列表。