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
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?
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
ok.... will work this ohter code out and post it here!....
THX
here we go
Link to code
Can anybody help me`?
,,mmm,,, dowesn't seem so!
....
We probably can butQuote:
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.
ok here's the code:
it should solve such a differential equation.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?
Also, you need to use more whitespace...Quote:
Originally posted by Mistert77
ok here's the code:
it should solve such a differential equation.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?
... 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.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
Your current formatting makes it so hard to read many of us aren't going to try.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.