Compute Jacobian matrix of function with multiple components

103 Views Asked by At

Hi I'm trying to fit a curve with gsl (Gnu Scientific Library). For the curve fit I need a Jacobian matrix, something I've never heard of before. I'm trying to wrap my head around it, but I simply don't understand how to translate the many tutorials on the topic to my curve.

A simplified version of the function is as follows:

$$ f(x, a, b, n) = \begin{cases} 0, & \text{for } x<a \newline (nx-na)/(b-a), & \text{for } a \leq x < b \newline n , & \text{for } b \leq x \newline \end{cases} $$

This function is 0 at the start, until $x$ is greater or equal to $a$, then it rises to $n$ and arrives there at $b$. If someone could explain to me how to compute the Jacobian matrix for that, it would already help a lot.

The full function is as follows: $$ f(x, a, b, c, d, n, m) = \begin{cases} 0, & \text{for } x<a \newline (nx-ca)/(b-a), & \text{for } a \leq x < b \newline (n - m)(x-b)/(c-b)+n , & \text{for } b \leq x < c \newline -m(x-c)/(d-c)+m, & \text{for } c \leq x < d \newline 0, & \text{for } d \leq x \end{cases} $$

Which would basically look like a trapezoid. See $n$ and $m$ as amplitudes of a "flat" top. and $a$ as the start, $b$ the start of the flat top, $c$ the end of the flat top, and $d$ as the end. For this I compute the following matrix with c. I'm not really sure how to translate that to math notation. I think I'm making a mistake somewhere

for (size_t i = 0; i < X; ++i) {
        double x = X[i];
        for(size_t j = 0; j < 6; j++) {
            gsl_matrix_set(J, i, j, 0);
        }
        if (x < a) {
            // do nothing, is already initialized to zero
        } else if (x < b) {
             gsl_matrix_set(J, i, 0, n*(x-b)/(pow(b-a,2)));  
             gsl_matrix_set(J, i, 1, n*(a-x)/(pow(b-a,2))); 
             gsl_matrix_set(J, i, 4, (x-start)/(b- a)); 
        } else if (x< c) {
            gsl_matrix_set(J, i, 1, (c-x)*(n-m)/(pow(c-b, 2))); 
            gsl_matrix_set(J, i, 2, (a-x)*(n-m)/(pow(c-b, 2))); 
            gsl_matrix_set(J, i, 4, (c-x)/(c- a)); 
            gsl_matrix_set(J, i, 5, (x-start)/(c- a)); 
        } else if (x < d) {
            gsl_matrix_set(J, i, 2, m*(d- x)/(pow(end-c,2))); 
            gsl_matrix_set(J, i, 3, m*(x-c)/(pow(d-c,2))); 
            gsl_matrix_set(J, i, 5, (d-x)/(d-c));
        } else {
            // do nothing, is already initialized to zero
        }
    }
}

```