C array of functions

Question

I have a problem with a series of functions. I have an array of 'return values' (i compute them through matrices) from a single function sys which depends on a integer variable, lets say, j, and I want to return them according to this j , i mean, if i want the equation number j, for example, i just write sys(j) For this, i used a for loop but i don't know if it's well defined, because when i run my code, i don't get the right values. Is there a better way to have an array of functions and call them in a easy way? That would make easier to work with a function in a Runge Kutta method to solve a diff equation.

I let this part of the code here: (c is just the j integer i used to explain before)

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
int N=3;
double s=10.;
//float r=28.;
double b=8.0/3.0;

/ * Define functions * /
double sys(int c,double r,double y[])
{
    int l,m,n,p=0;
    double tmp;
    double t[3][3]={0};
    double j[3][3]={{-s,s,0},{r-y[2],-1,-y[0]},{y[1],y[0],-b}}; //Jacobiano
    double id[3][3] = { {y[3],y[6],y[9]} , {y[4],y[7],y[10]} , {y[5],y[8],y[11]} };
    double flat[N*(N+1)];

    // Multiplication of matrices J * Y
    for(l=0;l<N;l++)
    {
        for(m=0;m<N;m++)
        {
            for(n=0;n<N;n++)
            {
                t[l][m] += j[l][n] * id[n][m];
            }
        }
    }

    // Transpose the matrix (J * Y) -> () t
    for(l=0;l<N;l++)
    {
        for(m=l+1;m<N;m++)
        {
            tmp = t[l][m];
            t[l][m] = t[m][l];
            t[m][l] = tmp;
        }
    }

    // We flatten the array to be left in one array 
    for(l=0;l<N;l++)
    {
        for(m=0;m<N;m++)
        {
            flat[p+N] = t[l][m];
        }
    }

    flat[0] = s*(y[1]-y[0]);
    flat[1] = y[0]*(r-y[2])-y[1];
    flat[2] = y[0]*y[1]-b*y[2];

    for(l=0;l<(N*(N+1));l++)
    {
        if(c==l)
        {
            return flat[c];
        }
    }
}

EDIT ----------------------------------------------------------------

Ok, this is the part of the code where i use the function

int main(){
output = fopen("lyapcoef.dat","w");
int j,k;
int N2 = N*N;
int NN = N*(N+1);
double r;
double rmax = 29;
double t = 0;
double dt = 0.05;
double tf = 50;
double z[NN]; // Temporary matrix for RK4
double k1[N2],k2[N2],k3[N2],k4[N2];
double y[NN]; // Matrix for all variables

/* Initial conditions */

double u[N];
double phi[N][N];
double phiu[N];
double norm;
double lyap;
//Here we integrate the system using Runge-Kutta of fourth order

for(r=28;r<rmax;r++){   
    y[0]=19;
    y[1]=20;
    y[2]=50;
    for(j=N;j<NN;j++) y[j]=0;
    for(j=N;j<NN;j=j+3) y[j]=1; // Identity matrix for y from 3 to 11 
    while(t<tf){
        /* RK4 step 1 */
        for(j=0;j<NN;j++){
            k1[j] = sys(j,r,y)*dt;
            z[j] = y[j] + k1[j]*0.5;
        }
        /* RK4 step 2 */
        for(j=0;j<NN;j++){
            k2[j] = sys(j,r,z)*dt;
            z[j] = y[j] + k2[j]*0.5;
        }
        /* RK4 step 3 */
        for(j=0;j<NN;j++){
            k3[j] = sys(j,r,z)*dt;
            z[j] = y[j] + k3[j];
        }
        /* RK4 step 4 */
        for(j=0;j<NN;j++){
            k4[j] = sys(j,r,z)*dt;
        }
        /* Updating y matrix with new values */
        for(j=0;j<NN;j++){
            y[j] += (k1[j]/6.0 + k2[j]/3.0 + k3[j]/3.0 + k4[j]/6.0);
        }
        printf("%lf %lf %lf \n",y[0],y[1],y[2]);
        t += dt; 
    }

Please translate your comments to english before posting next time, as well as improving your code formatting. — Richard J. Ross III
– Richard J. Ross III, Commented Mar 4, 2012 at 23:26
Did you really use the letter 'l' as a loop variable? Thats nearly impossible to read or support. I would suggest ignoring the mathematical part of this problem and write a simple use-case/example of the scenario and either debug it yourself or present your question without the complexity of the real world example. There are too many unknowns about what you are asking, and the other responders are latching on to responses which may not answer the actual question. — Paul Fox
– Paul Fox, Commented May 25, 2013 at 22:45

Gareth McCaughan · Accepted Answer · 2012-03-04 23:09:06Z

1

Since you're actually computing all these values at the same time, what you really want is for the function to return them all together. The easiest way to do this is to pass in a pointer to an array, into which the function will write the values. Or perhaps two arrays; it looks to me as if the output of your function is (conceptually) a 3x3 matrix together with a length-3 vector.

So the declaration of sys would look something like this:

void sys(double v[3], double JYt[3][3], double r, const double y[12]);

where v would end up containing the first three elements of your flat and JYt would contain the rest. (More informative names are probably possible.)

Incidentally, the for loop at the end of your code is exactly equivalent to just saying return flat[c]; except that if c happens not to be >=0 and <N*(N+1) then control will just fall off the end of your function, which in practice means that it will return some random number that almost certainly isn't what you want.

answered Mar 4, 2012 at 23:09

Gareth McCaughan

20k1 gold badge43 silver badges64 bronze badges

Sign up to request clarification or add additional context in comments.

5 Comments

PerroNoob Over a year ago

Mm yes, but i do need the values in a single array, that's why I "flatten" the matrix in one of the last steps in the sys function. The first three flat are the Lorenz system, the rest are anoher diff equations that depends on the first, so that's why i need to compute all the matrices again each time i call the function

PerroNoob Over a year ago

I was asking if it's well defined or if there's some error making a for loop with the return inside it. Or if there's a better way to get that same output, i just don't wanna write the 12 different equations in sepparated functions because that would make a huge runge Kutta in the main part of the code

Gareth McCaughan Over a year ago

If you want the outputs in a single flattened array, that's fine. (So: void sys(double v[12], double r, const double y[12]);.) Unless I've badly misunderstood the calling code you posted, you really truly don't need to keep calling sys repeatedly. For instance, in your "RK4 step 1" block, sys will compute the same stuff in flat each time. Make it print out the comments of flat and check if you don't believe me.

Gareth McCaughan Over a year ago

The for loop with a return inside it is pointless and inefficient and will make the function exhibit undefined behaviour if the value of c is out of range -- but, apart from that, there's nothing wrong with it :-). It's not illegal C or anything.

PerroNoob Over a year ago

Ok, so i call the function before each step of RK4 and then just use flat[] in the loop? That would be much faster . Thanks, i'll try it out.

Jonathan Leffler · Accepted Answer · 2012-03-05 02:59:39Z

1

Your function sys() does an O(N³) calculation to multiply two matrices, then does a couple of O(N²) operations, and finally selects a single number to return. Then it is called the next time and goes through most of the same processing. It feels a tad wasteful unless (even if?) the matrices are really small.

The final loop in the function is a little odd, too:

for(l=0;l<(N*(N+1));l++)
{
    if(c==l)
    {
        return flat[c];
    }
}

Isn't that more simply written as:

return flat[c];

Or, perhaps:

if (c < N * (N+1))
    return flat[c];
else
    ...do something on disastrous error other than fall off the end of the
    ...function without returning a value as the code currently does...

answered Mar 5, 2012 at 2:59

Jonathan Leffler

759k145 gold badges961 silver badges1.3k bronze badges

Comments

alexis · Accepted Answer · 2012-03-04 23:10:59Z

0

I don't see where you are selecting an algorithm by the value of j. If that's what you're trying to describe, in C you can have an array of pointers to functions; you could use a numerical index to choose a function from the array, but you can also pass a pointer-to-a-function to another function that will call it.

That said: Judging from your code, you should keep it simple. If you want to use a number to control which code gets executed, just use an if or switch statement.

switch (c) {
case 0:
  /* Algorithm 0 */
  break;
case 1:
  /* Algorithm 1 */
  etc.

edited Mar 4, 2012 at 23:10

answered Mar 4, 2012 at 23:01

alexis

50.4k18 gold badges107 silver badges172 bronze badges

4 Comments

ruakh Over a year ago

Re: "I don't see where you are selecting an algorithm by the value of j": it's at the very bottom, at the return statement. (The current code computes values using all three algorithms, but attempts to only return the value that was computed by the specified algorithm.)

alexis Over a year ago

Use a switch statement. For starters, you can replace the whole last loop with return flat[c].

ruakh Over a year ago

It sounds like you've confused me with the OP?

PerroNoob Over a year ago

Thanks for yours answers. If i use switch i would have to write like 14 case . I wanted to avoid that, that's the reason why i used a for loop . The matrices should always be computed, , the important part is at the bottom, like you said, so if i call the function using j, it throws me the j value of the array flat [] . I was wondering if it was well defined

Collectives™ on Stack Overflow

C array of functions

3 Answers 3

5 Comments

Comments

4 Comments

Your Answer

Hot Network Questions

Collectives™ on Stack Overflow

3 Answers 3

5 Comments

Comments

4 Comments

Your Answer

Sign up or log in

Post as a guest

Related