问题描述
给出
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 应该有一个合适的列表。