Hi,
I hope this kind of problem has already been faced by users of this
mail-list.
I have a code that solves a relatively simple equation, it is somewhat an
Helmotz equation (with negative parameter) whose solution is combination of
an exponential increasing+ exponential decreasing function
(one can treat first this very simply analytically).
This is the equation
y''[r]+2/r y'[r]-k y[r]=0
analytical solutions are
y[r]=C1 Exp[- sqrt(k) r] /r+ C2 Exp[ sqrt(k) r] /(Sqrt[k]* r)
If one is interested in only the exponential decreasing solution
can just simply set C2=0.
For a numerical solution, one must instead put initial conditions on the
function and its derivative obtained from the analytical treatment
(which is equivalent to dropping the exponential-increasing solution), so
that the numerical code should of course integrate and obtain only the
exponential decreasing solution (C2=0).
After doing that, depending on the method and accuracy, you get the
desired exponential decay but after a while, (for higher initial K-values)
at some point the exponential increasing solution starts to dominate.
It is of course only matter of numerical integration and not of wrong
initial conditions.
What I'm wondering is whether there is some way to force the code anyway to
avoid the fake exponential increasing, or it must be held by the user
(e.g., stopping when derivative becomes positive).
Presently I use
gsl_odeiv_control *k = gsl_odeiv_control_y_new(eps,eps);
gsl_odeiv_control *g = gsl_odeiv_control_yp_new(eps,eps);
For interested users, I also attach the code.
The command line after standard compilation is very simple, e.g.,
./test_hmz yi=1 kpar=0.2 meth=bsimp xmax=100 eps=1e-14
(for meth you can set bsimp, rk4, gear1, gear2, rk8pd, rkck),
yi is th initial value of the function (doesn't modify the general solution
behavior)
kpar is the parameter of the function (the higher kpar, the steepest the
exponential decay)
xmax is the right bound of the integration domain (which starts from x=1)
eps is the precision of the GSL controlling function.
You'll see the output written in the two-column file "results_hmz.dat",
whose plot immediately shows the behaviour.
I hope there are some tricks I don't know to get the right behavior.
Thank you in advance
Ruben
#include<stdlib.h>
#include<stdio.h>
#include <string.h>
#include <stdarg.h>
#include<math.h>
#include<gsl/gsl_errno.h>
#include<gsl/gsl_odeiv.h>
#include<gsl/gsl_const_cgs.h>
#include<gsl/gsl_integration.h>
#include<gsl/gsl_matrix.h>
double numpar (char * parameter);
char * strrev(char * string);
char* stringpar (char * parameter);
/* ##################################################################################### */
/* Solution of the Poisson equation with screened potential */
/* using spherical coordinates */
/* y''[r]+2/r y[r]- k y[r]=0 */
/* The linearly independent solutions are */
/* y[r]=C1 Exp[- sqrt(k) r]/r+ C2 Exp[ sqrt(k) r]/(Sqrt[k]* r) */
/* thus there is an exponential decreasing and exponential increasing function. */
/* If one is interested only the the first one, then C2=0 */
/* Solving the equation with initial conditions requires to */
/* set the function at the value r0, y[r0] and then its derivative */
/* at the value derived analytically by derivation of the solution with C2=0 */
/* #################################################################################### */
#define Pi 3.14159
/* Jacobian */
int jac (double x, const double y[], double *dfdy, double dfdt[], void *params) {
double kpar=((double *)params)[0];
double dyi=((double *)params)[1];
double x0=((double *)params)[2];
gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy,2,2);
gsl_matrix *m=&dfdy_mat.matrix;
gsl_matrix_set(m,0,0,0);
gsl_matrix_set(m,0,1,1);
gsl_matrix_set(m,1,0,kpar);
gsl_matrix_set(m,1,1,-2/x);
dfdt[0]=0;
dfdt[1]=2/pow(x0,2)*dyi;
printf("Jacobian y[0] %5.3e dfdt[0] %5.3e dfdt[1] %5.3e\n", y[0],dfdt[0], dyi);
return GSL_SUCCESS;
}
double xmax;
char *output;
int funct (double x, const double y[], double f[], void *params);
int main (int argc, char *argv[]) {
FILE *dat;
char *meth;
double x,eps, xmin, kpar, yi, dyi, c1, zio;
/* ######################################################################*/
if (argc < 2) {
printf("Test the solution of the Helmotz equation in spherical coordinates\n");
printf("test_algorithm yi=yi kpar=kpar meth=meth eps=eps [xmax=xmax]\n\n");
exit(1);
}
/* ######################################################################*/
int t;
for (t=1; t<argc; t++) {
if (strstr(argv[t], "kpar=") != NULL) kpar=numpar(argv[t]);
if (strstr(argv[t], "meth=") != NULL) meth=stringpar(argv[t]);
if (strstr(argv[t], "eps=") != NULL) eps=numpar(argv[t]);
if (strstr(argv[t], "xmax=") != NULL) xmax=numpar(argv[t]);
if (strstr(argv[t], "yi=") != NULL) yi=numpar(argv[t]);
/* if (strstr(argv[t], "dyi=") != NULL) dyi=numpar(argv[t]); */
if (strstr(argv[t], "zio=") != NULL) zio=numpar(argv[t]);
}
if (!kpar) {
printf("Parameter kpar not set\n");
return(1);
}
if (! yi) {
printf("Set the initial values for the function\n");
return(1);
}
if (!output) output="results_hmz.dat";
if (!eps) eps=1e-7;
if (!xmax) xmax=100.;
if (!meth) {
printf("Parameter meth not set\n");
return(1);
}
printf("kpar %5.3f meth %s eps %5.3e\n", kpar, meth, eps);
/* ######################################################################*/
/* Choose the integration methd */
/* ######################################################################*/
const gsl_odeiv_step_type *T;
if (strcmp(meth,"rk4") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rk4; */
T = gsl_odeiv_step_rk4;
} else if (strcmp(meth,"gear1") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_gear1; */
T = gsl_odeiv_step_gear1;
} else if (strcmp(meth,"gear2") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_gear2; */
T = gsl_odeiv_step_gear2;
} else if (strcmp(meth,"rkck") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_rkck;
} else if (strcmp(meth,"rk8pd") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_rk8pd;
} else if (strcmp(meth,"bsimp") ==0) {
/* const gsl_odeiv_step_type *T = gsl_odeiv_step_rkck; */
T = gsl_odeiv_step_bsimp;
} else {
printf("Allowed methods are: rk4, gear1, gear2, rkck, rk8pd, bsimp\n\n");
exit(1);
}
/* void *jac; */
gsl_odeiv_step *s = gsl_odeiv_step_alloc (T,2);
gsl_odeiv_control *k = gsl_odeiv_control_y_new(eps,eps);
gsl_odeiv_control *g = gsl_odeiv_control_yp_new(eps,eps);
gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (2);
/* printf("Receveing r0=%5.4e pc=%5.4e alpha=%5.4e gamma_val=%5.4e\n", r0, pc, alpha, gamma_val); */
/* ######################################################################## */
/* Set the initial value of the independent variable and stepsize */
double xinit=pow(10,-8);
xmin=1;
x=xmin;
double h=1e-4;
/* Determine the value of the constant for y[xmin]=yi */
c1=exp(sqrt(kpar)*xmin)*xmin*yi;
/* Compute value of the derivative */
dyi=-c1* exp(-sqrt(kpar)*xmin)/pow(xmin,2)- c1* exp(-sqrt(kpar)*xmin)*sqrt(kpar)/xmin;
dyi=-yi*(sqrt(kpar)+1/xmin);
/* ######################################################################## */
double params[3];
params[0]=kpar;
params[1]=dyi;
params[2]=xmin;
gsl_odeiv_system sys={funct, jac, 2, params};
printf("!Value c1 %5.3f yi %5.3f dyi %5.3e\n", c1,yi, dyi);
double y[2]={yi, dyi};
/* ######################################################################## */
/* Fixed step size integration */
/* double dydx_in[2], dydx_out[2], y_err[2]; */
/* GSL_ODEIV_FN_EVAL (&sys, t,y,dydx_in); */
/* ######################################################################## */
int status, count=0, count_loop=0;
double xprev=-10;
double minR=1000;
int countns=0;
dat=fopen(output, "w");
printf("!Value c1 %5.3f yi %5.3f dyi %5.3e x %5.3e xmax %5.3f\n", c1,y[0], y[1], x, xmax);
while (x < xmax) {
status = gsl_odeiv_evolve_apply (e,g,s,&sys, &x,xmax , &h, y);
/* status = gsl_odeiv_step_apply (s,x,h,y, y_err, dydx_in, dydx_out, &sys); */
if (x > xprev) {
fprintf(dat, "%5.4e %5.4e\n", x, y[0]);
}
if (status !=GSL_SUCCESS) {
break;
}
/* ############################################################################################ */
xprev=x;
count++;
/* dydx_in[0]=dydx_out[0]; */
/* dydx_in[1]=dydx_out[1]; */
/* x=x+h; */
}
fclose(dat);
gsl_odeiv_evolve_free(e);
gsl_odeiv_control_free(g);
gsl_odeiv_step_free(s);
}
/* ######################################################################## */
/*Field equations*/
/* ######################################################################## */
int funct (double x, const double y[], double f[], void *params) {
int static iter=0, iterns=0;
double kpar=((double *)params)[0];
f[0] = y[1];
f[1] =-2 * y[1]/x + kpar* y[0];
printf("Function x %5.3e y[0] %5.3e f[0] %5.3e\n", x, y[0], f[0]);
return GSL_SUCCESS;
}
/* ################################################################################################ */
char * strrev(char * string) {
int length = strlen(string);
char * result = malloc(length+1);
if( result != NULL ) {
int i,j; result[length] = '\0';
for ( i = length-1, j=0; i >= 0; i--, j++ ) result[j] = string[i];
}
return result;
}
/* ################################################################################################ */
char* stringpar (char * parameter) {
int t, ch='=';
double value;
char *file;
int substr=strlen(strrev(parameter))-strlen(strchr(strrev(parameter),ch));
char *dest=malloc(substr+1);
file=strrev((strncpy(dest, strrev(parameter), substr)));
return(file);
}
double numpar (char * parameter) {
int t, ch='=';
double value;
/* printf("%s\n", parameter); */
int substr=strlen(strrev(parameter))-strlen(strchr(strrev(parameter),ch));
char *dest=malloc(substr);
value=atof(strrev((strncpy(dest, strrev(parameter), substr))));
/* printf("%s m=%lf\n", strrev(strncpy(dest, strrev(parameter), substr)), m); */
return(value);
}
