Бисекция для корня функции (дигамма и вектор); C программирование - PullRequest
Новая
       13

Бисекция для корня функции (дигамма и вектор); C программирование

0 голосов
/ 20 марта 2019

Я хочу написать код ac, в котором он может найти корень (альфа) функции "f" на основе метода Bisection, заданного вектором X = [x_1, x_2, ..., x_n].Вектор и «n» (размер вектора) являются произвольными.Интересующая функция:

$$ f (\ alpha) = -n \ frac {\ Gamma '(\ alpha)} {\ Gamma (\ alpha)} + \ sum_ {i = 1} ^ {n} \ log (x_i). $$

Я не смог сделать вышеизложенное в формате формулы MathJax!Мой C-код выглядит следующим образом, но он не работает должным образом ... Любые решения приветствуются. 

#include<stdio.h>
#include<math.h>
#include <assert.h>

/* Digamma Function <Derivative of Log-Gamma> Based on the Series Expansion */
double digamma(double x){
    double result = 0, xx, xx2, xx4;
    assert(x > 0);
    for ( ; x < 7; ++x)
        result -= 1/x;
    x -= 1.0/2.0;
    xx = 1.0/x;
    xx2 = xx*xx;
    xx4 = xx2*xx2;
    result += log(x)+(1./24.)*xx2-(7.0/960.0)*xx4+(31.0/8064.0)*xx4*xx2-(127.0/30720.0)*xx4*xx4;
    return result;
}


float fun(float alpha){
    int n;  
    float *xvec, eq;
    n = sizeof(xvec);
    eq = -n*digamma(alpha)+sum(log(xvec));
    return eq;
}

void bisection (float *xvec, int n)
/* this function performs and prints the result of one iteration */
{
    float a = fmin(*xvec);
    float b = fmax(*xvec);
    float *x=(a+b)/2;
    int *itr;
    ++(*itr);
    printf("Iteration no. %3d X = %7.5f\n", *itr, *x);
}
void main ()
{
    int itr = 0, maxmitr = 200; /* Maximum Iteration */
    float x, x1;
    float allerr = 0.0001 /* Allowed Error */
    bisection (&xvec, &n);
    do
    {
        if (fun(a) * fun(x) < 0)
            b = x;
        else
            a = x;
        bisection (&x1, &n);
        if (fabs(x1-x) < allerr)
        {
            printf("After %d iterations, root = %6.4f\n", itr, x1);
            return 0;
        }
        x = x1;
    }
    while (itr < maxmitr);
    printf("The Solution Does not Converge!");
    return 0;
}
...