Stiff sets of equations

[Mistert77]

kennt sich damit jemand aus?

kennt das jemand?

ich müßte da was programmieren um so eine DGL zu lösen...und bräuchte Hilfe!

Hab hier das Buch "numerical Recipes in C" da sthet eigentlich alles drin wies geht und mit dem code... aber ich raffs nicht!

ich müßte das so programmieren, dass man eine solche DGL testen kann!

[Mistert77]

hier ist der 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");
}

weiß nicht wo die Gleichung reingeschrieben werden soll.