Thread: Stiff sets of equations

  1. #1
    Registered User
    Join Date
    Oct 2003
    Posts
    5

    Stiff sets of equations

    does anybody know about this?

    I have to make a program that can solve such a DGL!

    I got this book "numerical recipes in C", actually there is some kind of solution in it... but I do not really understand it.

    Is there anybody who can help me?

  2. #2
    Comment your source code! Lynux-Penguin's Avatar
    Join Date
    Apr 2002
    Posts
    533
    explain the problem more accurately. Write some code or psuedo code on how YOU think it should be.

    You need to show us some effort or some ideas to work on, we really despise doing other people's homework, especially when we're here instead of doing our own =/

    -LC
    Asking the right question is sometimes more important than knowing the answer.
    Please read the FAQ
    C Reference Card (A MUST!)
    Pointers and Memory
    The Essentials
    CString lib

  3. #3
    Registered User
    Join Date
    Oct 2003
    Posts
    5
    ok.... will work this ohter code out and post it here!....
    THX

  4. #4
    Registered User
    Join Date
    Oct 2003
    Posts
    5
    here we go

    Link to code

  5. #5
    Registered User
    Join Date
    Oct 2003
    Posts
    5
    Can anybody help me`?

    ,,mmm,,, dowesn't seem so!


    ....

  6. #6
    Been here, done that.
    Join Date
    May 2003
    Posts
    1,164
    Originally posted by Mistert77
    Can anybody help me`?

    ,,mmm,,, dowesn't seem so!


    ....
    We probably can but
    1) most of us have lives outside of this board so we're not sitting here waiting for you to post your question. When we log in, we do what we can, and that may take a few hours.
    2) When you post code, the code must be visible, not hidden elsewhere on the web. Post code means "Post the Code" not give a link.
    3) You really need to specify what your problem is. Cryptic messages like "Does anyone know about this?" are not helpful to our understanding what problem you are having so we cannot help -- yet.

    Please post code and explain what is and is not working properly.
    Definition: Politics -- Latin, from
    poly meaning many and
    tics meaning blood sucking parasites
    -- Tom Smothers

  7. #7
    Registered User
    Join Date
    Oct 2003
    Posts
    5
    ok here's the code:

    Code:
    #include <math.h>
    #include "nrutil.h"
    #define SAFETY 0.9
    #define GROW 1.5
    #define PGROW -0.25
    #define SHRNK 0.5
    #define PSHRNK (-1.0/3.0)
    #define ERRCON 0.1296
    #define MAXTRY 40
    #define GAM (1.0/2.0)
    #define A21 2.0
    #define A31 (48.0/25.0)
    #define A32 (6.0/25.0)
    #define C21 -8.0
    #define C31 (372.0/25.0)
    #define C32 (12.0/5.0)
    #define C41 (-112.0/125.0)
    #define C42 (-54.0/125.0)
    #define C43 (-2.0/5.0)
    #define B1 (19.0/9.0)
    #define B2 (1.0/2.0)
    #define B3 (25.0/108.0)
    #define B4 (125.0/108.0)
    #define E1 (17.0/54.0)
    #define E2 (7.0/36.0)
    #define E3 0.0
    #define E4 (125.0/108.0)
    #define C1X (1.0/2.0)
    #define C2X (-3.0/2.0)
    #define C3X (121.0/50.0)
    #define C4X (29.0/250.0)
    #define A2X 1.0
    #define A3X (3.0/5.0)
    void stiff(float y[], float dydx[], int n, float *x, float htry, float eps,
    float yscal[], float *hdid, float *hnext,
    void (*derivs)(float, float [], float []))
    {
    void jacobn(float x, float y[], float dfdx[], float **dfdy, int n);
    void lubksb(float **a, int n, int *indx, float b[]);
    void ludcmp(float **a, int n, int *indx, float *d);
    int i,j,jtry,*indx;
    float d,errmax,h,xsav,**a,*dfdx,**dfdy,*dysav,*err;
    float *g1,*g2,*g3,*g4,*ysav;
    indx=ivector(1,n);
    a=matrix(1,n,1,n);
    dfdx=vector(1,n);
    dfdy=matrix(1,n,1,n);
    dysav=vector(1,n);
    err=vector(1,n);
    g1=vector(1,n);
    g2=vector(1,n);
    g3=vector(1,n);
    g4=vector(1,n);
    ysav=vector(1,n);
    xsav=(*x); 
    for (i=1;i<=n;i++) {
    ysav[i]=y[i];
    dysav[i]=dydx[i];
    }
    jacobn(xsav,ysav,dfdx,dfdy,n);
    h=htry;
    for (jtry=1;jtry<=MAXTRY;jtry++) {
    for (i=1;i<=n;i++) { 
    for (j=1;j<=n;j++) a[i][j] = -dfdy[i][j];
    a[i][i] += 1.0/(GAM*h);
    }
    ludcmp(a,n,indx,&d); 
    for (i=1;i<=n;i++) 
    g1[i]=dysav[i]+h*C1X*dfdx[i];
    lubksb(a,n,indx,g1); 
    for (i=1;i<=n;i++) 
    y[i]=ysav[i]+A21*g1[i];
    *x=xsav+A2X*h;
    (*derivs)(*x,y,dydx); 
    for (i=1;i<=n;i++) 
    g2[i]=dydx[i]+h*C2X*dfdx[i]+C21*g1[i]/h;
    lubksb(a,n,indx,g2); 
    for (i=1;i<=n;i++) 
    y[i]=ysav[i]+A31*g1[i]+A32*g2[i];
    
    *x=xsav+A3X*h;
    (*derivs)(*x,y,dydx); 
    for (i=1;i<=n;i++) 
    g3[i]=dydx[i]+h*C3X*dfdx[i]+(C31*g1[i]+C32*g2[i])/h;
    lubksb(a,n,indx,g3); 
    for (i=1;i<=n;i++) 
    g4[i]=dydx[i]+h*C4X*dfdx[i]+(C41*g1[i]+C42*g2[i]+C43*g3[i])/h;
    lubksb(a,n,indx,g4); 
    for (i=1;i<=n;i++) { 
    y[i]=ysav[i]+B1*g1[i]+B2*g2[i]+B3*g3[i]+B4*g4[i];
    err[i]=E1*g1[i]+E2*g2[i]+E3*g3[i]+E4*g4[i];
    }
    *x=xsav+h;
    if (*x == xsav) nrerror("stepsize not significant in stiff");
    errmax=0.0; 
    for (i=1;i<=n;i++) errmax=FMAX(errmax,fabs(err[i]/yscal[i]));
    errmax /= eps; 
    if (errmax <= 1.0) { 
    *hdid=h;
    *hnext=(errmax > ERRCON ? SAFETY*h*pow(errmax,PGROW) : GROW*h);
    free_vector(ysav,1,n);
    free_vector(g4,1,n);
    free_vector(g3,1,n);
    free_vector(g2,1,n);
    free_vector(g1,1,n);
    free_vector(err,1,n);
    free_vector(dysav,1,n);
    free_matrix(dfdy,1,n,1,n);
    free_vector(dfdx,1,n);
    free_matrix(a,1,n,1,n);
    free_ivector(indx,1,n);
    return;
    
    } else { 
    *hnext=SAFETY*h*pow(errmax,PSHRNK);
    h=(h >= 0.0 ? FMAX(*hnext,SHRNK*h) : FMIN(*hnext,SHRNK*h));
    }
    } 
    nrerror("exceeded MAXTRY in stiff");
    }
    it should solve such a differential equation.

    I don't know here to put this equation?

  8. #8
    Been here, done that.
    Join Date
    May 2003
    Posts
    1,164
    Originally posted by Mistert77
    ok here's the code:

    Code:
    .... snip ....
    #define C3X (121.0/50.0)
    #define C4X (29.0/250.0)
    #define A2X 1.0
    #define A3X (3.0/5.0)
    void stiff(float y[], float dydx[], int n, float *x, float htry, float eps,
    float yscal[], float *hdid, float *hnext,
    void (*derivs)(float, float [], float []))
    {
    // These statements look like prototypes:
    void jacobn(float x, float y[], float dfdx[], float **dfdy, int n);
    void lubksb(float **a, int n, int *indx, float b[]);
    void ludcmp(float **a, int n, int *indx, float *d);
    // They have to be placed outside of all function, below your defines
    
    int i,j,jtry,*indx;
    float d,errmax,h,xsav,**a,*dfdx,**dfdy,*dysav,*err;
    float *g1,*g2,*g3,*g4,*ysav;
    indx=ivector(1,n);
    a=matrix(1,n,1,n);
    dfdx=vector(1,n);
    dfdy=matrix(1,n,1,n);
    dysav=vector(1,n);
    err=vector(1,n);
    g1=vector(1,n);
    g2=vector(1,n);
    g3=vector(1,n);
    .... more snipping ....
    it should solve such a differential equation.

    I don't know here to put this equation?
    Also, you need to use more whitespace...
    Code:
    for (i = 1; i <= n; i++) 
    g3[i]  = dydx[i] + h * C3X * dfdx[i] + (C31 * g1[i] + C32 * g2[i]) / h;
    lubksb(a, n, indx, g3);
    ... and use more parentheses to combine the proper elements. You have a mixture of mults & adds that if not combined properly will give you the wrong answer.

    For example,
    a = 3 * 5+4
    The result is 19 when the desired result is 27. The addition must be performed first so you need to tell the compiler:
    a = 3 * (5+4)

    And indent blocks of code
    Code:
    jacobn(xsav, ysav, dfdx, dfdy, n);
    h = htry;
    for (jtry = 1; jtry <= MAXTRY; jtry++) 
    {
        for (i = 1; i <= n; i++) 
        { 
            for (j=1;j<=n;j++) a[i][j] = -dfdy[i][j];
            a[i][i] += 1.0 / (GAM * h);
        }
        ludcmp(a,n,indx,&d); 
        for (i = 1; i <= n; i++) 
            g1[i] = dysav[i] + h * C1X * dfdx[i];
        lubksb(a, n, indx, g1); 
        for (i = 1;i <= n; i++) 
            y[i] = ysav[i] + A21 * g1[i];
        *x = xsav + A2X * h;
    Your current formatting makes it so hard to read many of us aren't going to try.

    As for where you put the equation, I don't understand the question. The program is supposed to be a solution to the equation. You probably want to put in a main() block to ask the user the values of equation parameters, then call the function that solves it.

    Many (most?) of us are programmers, not math genuises (I can't spell mathematician) and don't understand what let alone how to solve a differential equation. You have to help us understand to process for us to check for your errors.
    Definition: Politics -- Latin, from
    poly meaning many and
    tics meaning blood sucking parasites
    -- Tom Smothers

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Problem aligning floating point numbers
    By esbo in forum C Programming
    Replies: 4
    Last Post: 01-05-2009, 08:09 PM
  2. Major Problems with GetSaveFileName() and GetOpenFileName()
    By CodeHacker in forum Windows Programming
    Replies: 8
    Last Post: 07-12-2004, 11:05 AM
  3. creating new sets
    By axon in forum C++ Programming
    Replies: 7
    Last Post: 12-03-2003, 06:37 PM
  4. LUP Decomposition (simultaneous equations)
    By Markallen85 in forum C Programming
    Replies: 6
    Last Post: 08-24-2003, 02:08 AM
  5. Just one Question?
    By Irish-Slasher in forum C++ Programming
    Replies: 6
    Last Post: 02-12-2002, 10:19 AM