问题描述
首先,我是 C++ 的新手。 Tbh 我发现很难“习惯”,但由于计算时间要求,上周我一直在尝试将脚本从 Python“翻译”到 C++。
首先是来自 Python 的简单 MWE: 试图解决 ax^2+bx+c = 0 :
def f(x,a,b,c):
return a*x**2+b*x+c
optimize.fsolve(f,-1.0e2,args=(1,1,-1))
根据起始的猜测符号返回 -1.618
或 0.618
好的,我在网上搜索的 C++ 有点复杂,使用了 GSL Root finding 库。
对于简单的非多项式函数,它的作用就像一个魅力,但是当二阶出现时, 端点的添加在您简化搜索“快速”解决方案时似乎是一个问题:
#include <stdio.h>
#include <math.h>
#include <functional>
#include <stdlib.h>
#include <iostream>
#include <vector>
#include <gsl/gsl_math.h>
#include <gsl/gsl_interp2d.h>
#include <gsl/gsl_spline2d.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_spline.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_roots.h>
struct my_f_params { double a; double b; double c; };
double
my_f (double x,void * p)
{
struct my_f_params * params = (struct my_f_params *)p;
double a = (params->a);
double b = (params->b);
double c = (params->c);
return (a * x + b) * x + c;
}
double root (struct my_f_params prms,double r)
{
int status;
int iter = 0,max_iter = 50;
const gsl_root_fsolver_type *T;
gsl_root_fsolver *s;
double x_lo= -5e0,x_hi=11e0;
gsl_function F;
F.function = &my_f;
F.params = &prms;
T = gsl_root_fsolver_falsepos;
s = gsl_root_fsolver_alloc (T);
gsl_root_fsolver_set (s,&F,x_lo,x_hi);
do
{
iter++;
gsl_set_error_handler_off();
status = gsl_root_fsolver_iterate (s);
r = gsl_root_fsolver_root (s);
x_lo = gsl_root_fsolver_x_lower (s);
x_hi = gsl_root_fsolver_x_upper (s);
status = gsl_root_test_interval (x_lo,x_hi,0.001);
printf("%f %f\n",x_hi);
}
while (status == GSL_CONTINUE && iter < max_iter);
return r;
}
int main(int argc,char const *argv[])
{
struct my_f_params params = {1,-1};
printf("root of x2+x1-1=0 %f\n",root(params,1.25));
return 0;
}
现在,如果起始 x_lo,x_hi
覆盖了 2 个解积分,它不会继续寻找最接近的一个,给出错误
gsl: falsepos.c:74: ERROR: endpoints do not sTraddle y=0
Default GSL error handler invoked.
Aborted (core dumped)
在我在这里发布之前已经尝试了很多来自 google 的东西。
真的非常感谢您的时间,任何事情都非常感谢!
解决方法
您需要的是研究这个类似的问题,并很好地解释如何使其工作error in GSL - root finding。 为了提供帮助,您可以使用 https://www.desmos.com/calculator/zuaqvcvpbz 来查看该函数,并且您可以像这样设置初始值: 双 x_lo = 0.0,x_hi = 1.0; 在您的实现中使其运行; 添加此代码将有助于为 x_lo 和 x_hi 找到合适的值:
gsl_set_error_handler_off(); // this turns off error reporting
int check = gsl_root_fsolver_set(s,&F,x_lo,x_hi);
if (check == GSL_EINVAL) {// this is the error code you got
do {
x_lo += 0.1; // it would be appropriate to check the sign in both
x_hi -= 0.1; // cases,to make sure interval is adjusted properly
check = gsl_root_fsolver_set(s,x_hi);
} while (check != 0);
}