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?

Printable View

- 10-13-2003Mistert77Stiff 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? - 10-13-2003Lynux-Penguin
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 - 10-14-2003Mistert77
ok.... will work this ohter code out and post it here!....

THX - 10-14-2003Mistert77
here we go

Link to code - 10-16-2003Mistert77
Can anybody help me`?

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

.... - 10-17-2003WaltPQuote:

*Originally posted by Mistert77*

**Can anybody help me`?**

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

....

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. - 10-27-2003Mistert77
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");

}

I don't know here to put this equation? - 10-27-2003WaltPQuote:

*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 ....

I don't know here to put this equation?

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);

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;

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.