使用C++库GSL解非线性方程

解方程是一个比较常见的功能,matlab,python都有简单方法实现。C++也有不少实现方法,下面就介绍一下GSL解方程的方法,代码参考了GSL源码中的demo。
GSL解方程的功能相对matlab要简单一些,目前看只支持实数解,并且只有一个解。当然对于大部分实际应用 已经足够了

假设需要解的方程组是如下2元非线性方程(方程右边不为0,移项一下就好了)

解方程的代码如下:
1.加入头文件

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

2.定义一个struct 保存可变常数

struct rparams
{
	double a;
	double b;
};

3.定义方程

int rosenbrock_f(const gsl_vector * x, void *params,
	gsl_vector * f)
{
	double a = ((struct rparams *) params)->a; //常数a
	double b = ((struct rparams *) params)->b; //常数b

	const double x0 = gsl_vector_get(x, 0); //未知数x
	const double x1 = gsl_vector_get(x, 1); //未知数y

	const double y0 = a * (4 - x0) + x1*x1 - b; //方程1
	const double y1 = b * (x1 - sqrt(x0)); //方程2

	gsl_vector_set(f, 0, y0); //set 方程1
	gsl_vector_set(f, 1, y1); //set 方程2

	return GSL_SUCCESS;
}

4.定义过程打印输出

int print_state_f(size_t iter, gsl_multiroot_fsolver * s)
{
	printf("iter = %3u x = % .3f % .3f "
		"f(x) = % .3e % .3e\n",
		iter,
		gsl_vector_get(s->x, 0),
		gsl_vector_get(s->x, 1),
		gsl_vector_get(s->f, 0),
		gsl_vector_get(s->f, 1));
	return 0;
}

4.定义主函数
详细说明看注释

int main(void)
{
	const gsl_multiroot_fsolver_type *T; 
	gsl_multiroot_fsolver *s;

	int status;
	size_t i, iter = 0;

	const size_t n = 2; // n 定义了元的多少
	struct rparams p = { 8.52550, 10.5665 }; // 定义常数a b 如果需要多次计算不同a b时的方程解,改这个p就可以了
	gsl_multiroot_function f = { &rosenbrock_f,n, &p }; // 定义解方程的函数

	double x_init[2] = { 10.0, 3.0 }; // 定义方程解预设值,这个值必须在真实解的附近,不然会解不出来。这个也是GSL的局限,我暂时没有找到解决的办法
	gsl_vector *x = gsl_vector_alloc(n);

	gsl_vector_set(x, 0, x_init[0]);
	gsl_vector_set(x, 1, x_init[1]);

	T = gsl_multiroot_fsolver_hybrids;
	s = gsl_multiroot_fsolver_alloc(T, n);
	gsl_multiroot_fsolver_set(s, &f, x);

	print_state_f(iter, s);

	do //这个循环迭代解方程,最多迭代1000次
	{
		iter++;
		status = gsl_multiroot_fsolver_iterate(s);
		print_state_f(iter, s);
		if (status)
			break;
		status = gsl_multiroot_test_residual(s->f, 1e-7);//判断解是否是真实解
	} while (status == GSL_CONTINUE && iter < 1000);

	printf("status = %s\n", gsl_strerror(status));

	gsl_multiroot_fsolver_free(s);
	gsl_vector_free(x);
	return 0;
}

5.运行结果
可以看到运行7次后得到解。可以修改预设的解试试,当预设解偏离太多真实解时,就没办法得到解了。

iter =   0 x =  10.000  3.000 f(x) = -5.272e+01 -1.715e+00
iter =   1 x =  3.171  2.082 f(x) =  8.419e-01  3.189e+00
iter =   2 x =  3.039  1.750 f(x) =  6.935e-01  7.627e-02
iter =   3 x =  3.120  1.760 f(x) =  3.637e-02 -6.537e-02
iter =   4 x =  3.128  1.769 f(x) = -4.847e-03  4.300e-03
iter =   5 x =  3.127  1.768 f(x) = -6.137e-05  4.227e-05
iter =   6 x =  3.127  1.768 f(x) = -2.508e-07  2.030e-07
iter =   7 x =  3.127  1.768 f(x) =  7.090e-10 -5.579e-10
status = success

6.使用导数,加快运算速度
可以手动定义方程组的偏导数,来加快计算过程。代码如下:
大部分内容一样的,多了一步计算偏导数的过程,详见注释

#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>

struct rparams
{
	double a;
	double b;
};

int
rosenbrock_f(const gsl_vector * x, void *params,
	gsl_vector * f)
{
	double a = ((struct rparams *) params)->a;
	double b = ((struct rparams *) params)->b;

	const double x0 = gsl_vector_get(x, 0);
	const double x1 = gsl_vector_get(x, 1);

	const double y0 = a * (4 - x0) + x1*x1 - b;
	const double y1 = b * (x1 - sqrt(x0));

	gsl_vector_set(f, 0, y0);
	gsl_vector_set(f, 1, y1);

	return GSL_SUCCESS;
}

int rosenbrock_df(const gsl_vector * x, void *params,
	gsl_matrix * J) //定义偏导数
{
	const double a = ((struct rparams *) params)->a;
	const double b = ((struct rparams *) params)->b;

	const double x0 = gsl_vector_get(x, 0); //未知数x
	const double x1 = gsl_vector_get(x, 1); //未知数y

	const double df00 = -a; // 方程1 对x的偏导数
	const double df01 = 2 * x1; // 方程1 对y的偏导数
	const double df10 = - 0.5 * b  /  sqrt(x0); // 方程2对x的偏导数
	const double df11 = b; //方程2 对y的偏导数

	gsl_matrix_set(J, 0, 0, df00);
	gsl_matrix_set(J, 0, 1, df01);
	gsl_matrix_set(J, 1, 0, df10);
	gsl_matrix_set(J, 1, 1, df11);

	return GSL_SUCCESS;
}

int rosenbrock_fdf(const gsl_vector * x, void *params,
	gsl_vector * f, gsl_matrix * J)
{
	rosenbrock_f(x, params, f);
	rosenbrock_df(x, params, J);
	return GSL_SUCCESS;
}

int print_state_fdf(size_t iter, gsl_multiroot_fdfsolver * s)
{
	printf("iter = %3u x = % .3f % .3f "
		"f(x) = % .3e % .3e\n",
		iter,
		gsl_vector_get(s->x, 0),
		gsl_vector_get(s->x, 1),
		gsl_vector_get(s->f, 0),
		gsl_vector_get(s->f, 1));
	return 0;
}

int main(void)
{
	const gsl_multiroot_fdfsolver_type *T;
	gsl_multiroot_fdfsolver *s;

	int status;
	size_t i, iter = 0;

	const size_t n = 2;
	struct rparams p = { 8.52550, 10.5665 };
	gsl_multiroot_function_fdf f = {&rosenbrock_f,&rosenbrock_df,&rosenbrock_fdf,n, &p };

	double x_init[2] = { 10.0, 3.0 };
	gsl_vector *x = gsl_vector_alloc(n);

	gsl_vector_set(x, 0, x_init[0]);
	gsl_vector_set(x, 1, x_init[1]);

	T = gsl_multiroot_fdfsolver_gnewton;
	s = gsl_multiroot_fdfsolver_alloc(T, n);
	gsl_multiroot_fdfsolver_set(s, &f, x);

	print_state_fdf(iter, s);

	do
	{
		iter++;
		status = gsl_multiroot_fdfsolver_iterate(s);
		print_state_fdf(iter, s);
		if (status)
			break;
		status = gsl_multiroot_test_residual(s->f, 1e-7);
	} while (status == GSL_CONTINUE && iter < 1000);

	printf("status = %s\n", gsl_strerror(status));

	gsl_multiroot_fdfsolver_free(s);
	gsl_vector_free(x);
	return 0;
}

输出:
可以看到收敛的更快了一点

iter =   0 x =  10.000  3.000 f(x) = -5.272e+01 -1.715e+00
iter =   1 x =  3.171  2.082 f(x) =  8.419e-01  3.189e+00
iter =   2 x =  3.114  1.765 f(x) =  1.009e-01  7.523e-04
iter =   3 x =  3.127  1.768 f(x) =  1.384e-05  4.293e-05
iter =   4 x =  3.127  1.768 f(x) =  1.667e-11  2.346e-15
status = success

发表回复

您的电子邮箱地址不会被公开。