GSL 求解非线性方程

侃侃_天下

于 2023-07-30 08:35:58 发布

阅读量269

点赞数 1

分类专栏： C++ 文章标签：算法 c++

本文链接：https://blog.csdn.net/douzi949389/article/details/132003322

版权

C++ 专栏收录该内容

76 篇文章 0 订阅

订阅专栏

#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <vector>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_deriv.h>
struct rparams
{
    size_t n;
};

int rosenbrock_f(const gsl_vector* x, void* params, gsl_vector* f)
{
    const double n = ((struct rparams*)params)->n;
    for(size_t i=0; i<n; ++i)
    {
        double xi = gsl_vector_get(x, i);
        double yi = sin(xi)- 0.5;
        gsl_vector_set(f, i, yi);
    }
    return GSL_SUCCESS;
}

int main(void)
{
    const gsl_multiroot_fsolver_type* T;
    gsl_multiroot_fsolver* s;
    std::vector<double> x_init = {1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0};
    size_t n = x_init.size();
    struct rparams p = { n };
    gsl_multiroot_function f = { &rosenbrock_f, n, &p };
//    gsl_function_struct g;
    gsl_vector* x = gsl_vector_alloc(n);
    for(size_t i=0; i<x_init.size(); ++i)
    {
        gsl_vector_set(x, i, x_init[i]);
    }

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

    int status = GSL_CONTINUE;
    for(size_t i=0; (status == GSL_CONTINUE) && i < 1000; ++i)
    {
        status = gsl_multiroot_fsolver_iterate(s);
        if(status)
        {
            break;
        }
        status = gsl_multiroot_test_residual(s->f, 1e-7);
    }

    for(size_t i=0; i<n; ++i)
    {
        printf ("% 7.4f", gsl_vector_get (s->x, i));
    }
    printf("\nstatus = %s\n", gsl_strerror(status));
    gsl_multiroot_fsolver_free(s);
    gsl_vector_free(x);
    return 0;
}

// 0.5236    2.6180    2.6180    2.6180    6.8068    6.8068    6.8068    8.9012   8.9012
// 0.5236    2.6180    2.6180    2.6180   13.0900    6.8068    6.8068    6.8068   8.9012
// 0.5236    2.6180    2.6180    2.6180   13.0900    6.8068    6.8068    2.6180   8.9012