Hello,
I have written a code to solve 5 non-linear coupled equations using gsl. I
have a few constants that are set in the beginning of the code and one
input parameter, rho, that is changed in a loop. The program works fine for
a given set of parameters (G300) but for a different one (GM4) it is
working only for small rho. I always use the solution found for rho for the
initial guess of rho+Delta_rho but I start receiving the "not making
progress towards the solution" error from a certain value of rho when using
the GM4 set of parameters.
In main, I have
gsl_multiroot_function f = {&Walecka_f, n, &p};
double x_init[5] = {gss_i, gww_i, grr_i, mun_i, mue_i};
gsl_vector *x=gsl_vector_alloc(n);
gsl_vector_set(x, 0, x_init[0]);
gsl_vector_set(x, 1, x_init[1]);
gsl_vector_set(x, 2, x_init[2]);
gsl_vector_set(x, 3, x_init[3]);
gsl_vector_set(x, 4, x_init[4]);
T = gsl_multiroot_fsolver_hybrid;
s = gsl_multiroot_fsolver_alloc(T, 5);
gsl_multiroot_fsolver_set(s, &f, x);
do
{
iter++;
status = gsl_multiroot_fsolver_iterate (s);
if (status) /* check if solver is stuck */
break;
status = gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
I I have also attached the code so if anyone can help, I would really
appreciate it.
Thanks,
Laura
/****************************************************************************/
/* This code was implemented for calculating nuclear EoS based on the */
/* Walecka model with matter composed of n, p, e, muons, Lambda, Sigma */
/* and Cascade */
/****************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_vector.h>
#include <gsl/gsl_multiroots.h>
#define hc 197.327
#define nmax 250000
#define epsilon 1e-8
#define me 0.511/hc // electron mass
#define mmu 105.0/hc // muon mass
#define mp 938.3/hc // proton mass
#define mn 939.6/hc // neutron mass
#define ml 1116.0/hc // Lambda mass
#define msp 1189.0/hc // Sigma + mass
#define msm 1197.0/hc // Sigma - mass
#define ms0 1193.0/hc // Sigma 0 mass
#define mc0 1315.0/hc //Chi 0 mass
#define mcm 1321.0/hc //Chi - mass
struct rparams
{
double gmw2; // (g_omega/m_omega)^2
double gms2; // (g_sigma/m_sigma)^2
double gmr2; // (g_rho/m_rho)^2
double b;
double c;
double xs; // x_sigma
double xw; // x_omega
double xr; // x_rho
double rho; //density
};
int Walecka_f (const gsl_vector * x, void *params, gsl_vector * f);
double Eq1 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq2 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq3 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq4 (double x1, double x2, double x3, double x4, double x5, void *params);
double Eq5 (double x1, double x2, double x3, double x4, double x5, void *params);
double kf2(double mub, double mue, double m8, double Vw, double Vr,double qe, double I3, void *params);
double dens_s(double m8, double k);
double dens(double k);
int main (void)
{
int status;
double mue;
double gs_s, gww, grr, b, c;
double mub, Vw, Vr, Vs;
double rhoaux;
double gms2, gmw2, gmr2;
double gss_i, mun_i, mue_i, gww_i, grr_i;
// double xs=1.0, xr=1.0, xw=1.0; //for G300
double xs=0.9, xr=0.9, xw=0.9; //for GM4
size_t iter = 0;
const size_t n = 5;
const gsl_multiroot_fsolver_type *T;
gsl_multiroot_fsolver *s;
rhoaux=0.01; //fm^-3
/* gmw2=4.733; //for G300
gms2=9.031;
gmr2=4.825;
b=0.003305;
c=0.01529;
*/
gmw2=7.149; //for GM4
gms2=11.79;
gmr2=4.411;
b=0.002947;
c=-0.001070;
gss_i=sqrt(300.0/(0.9999*939.0));
gss_i=asin(gss_i);
mun_i=mn;
mue_i=0.12*938.0/hc*pow(rhoaux/0.17, 2.0/3.0);
gww_i= sqrt(226.0/hc);
grr_i=1000.0/hc*rhoaux;
for (rhoaux=0.01; rhoaux<=1.7;){
iter=0;
struct rparams p = {gmw2, gms2, gmr2, b, c, xs, xw, xr, rhoaux};
gsl_multiroot_function f = {&Walecka_f, n, &p};
double x_init[5] = {gss_i, gww_i, grr_i, mun_i, mue_i};
gsl_vector *x=gsl_vector_alloc(n);
gsl_vector_set(x, 0, x_init[0]);
gsl_vector_set(x, 1, x_init[1]);
gsl_vector_set(x, 2, x_init[2]);
gsl_vector_set(x, 3, x_init[3]);
gsl_vector_set(x, 4, x_init[4]);
T = gsl_multiroot_fsolver_hybrid;
s = gsl_multiroot_fsolver_alloc(T, 5);
gsl_multiroot_fsolver_set(s, &f, x);
do
{
iter++;
printf ("iter = %3zu x = %g %g %g %g %g " "f(x) = %g %g %g %g %g \n", iter,
gsl_vector_get (s->x, 0),
gsl_vector_get (s->x, 1),
gsl_vector_get (s->x, 2),
gsl_vector_get (s->x, 3),
gsl_vector_get (s->x, 4),
gsl_vector_get (s->f, 0),
gsl_vector_get (s->f, 1),
gsl_vector_get (s->f, 2),
gsl_vector_get (s->f, 3),
gsl_vector_get (s->f, 4));
// print_state (iter, s);
status = gsl_multiroot_fsolver_iterate (s);
if (status) /* check if solver is stuck */
break;
status = gsl_multiroot_test_residual (s->f, 1e-7);
}
while (status == GSL_CONTINUE && iter < 1000);
printf ("status = %s\n", gsl_strerror (status));
Vs=gsl_vector_get (s->x,0);
Vw=gsl_vector_get (s->x,1);
Vr=gsl_vector_get (s->x,2);
mub=gsl_vector_get (s->x,3);
mue=gsl_vector_get (s->x,4);
//Mapeio:
gs_s=0.9999*mn*sin(Vs)*sin(Vs);
gww=Vw*Vw;
grr=Vr;
gsl_multiroot_fsolver_free (s);
gsl_vector_free (x);
gss_i=sqrt(gs_s/(0.9999*mn));
gss_i=asin(gss_i);
mun_i=mub;
mue_i=mue;
gww_i= sqrt(gww);
grr_i=grr;
if(rhoaux<0.1)
rhoaux+=0.005;
else
rhoaux+=0.025;
}
return 0;
}
int Walecka_f (const gsl_vector * x, void *params, gsl_vector * f)
{
const double x0 = gsl_vector_get (x, 0);
const double x1 = gsl_vector_get (x, 1);
const double x2 = gsl_vector_get (x, 2);
const double x3 = gsl_vector_get (x, 3);
const double x4 = gsl_vector_get (x, 4);
const double y0 = Eq1(x0, x1, x2, x3, x4, params);
const double y1 = Eq2(x0, x1, x2, x3, x4, params);
const double y2 = Eq3(x0, x1, x2, x3, x4, params);
const double y3 = Eq4(x0, x1, x2, x3, x4, params);
const double y4 = Eq5(x0, x1, x2, x3, x4, params);
gsl_vector_set (f, 0, y0);
gsl_vector_set (f, 1, y1);
gsl_vector_set (f, 2, y2);
gsl_vector_set (f, 3, y3);
gsl_vector_set (f, 4, y4);
return GSL_SUCCESS;
}
double Eq1 (double x1, double x2, double x3, double x4, double x5, void *params) //equacao para o campo sigma
{
double gms2 = ((struct rparams *) params)->gms2;
double b = ((struct rparams *) params)->b;
double c = ((struct rparams *) params)->c;
double xs = ((struct rparams *) params)->xs;
double Vs, Vr, Vw, mub, mue;
double rho_s=0.0, m8, k2, k;
double f;
//Mapeio:
Vs=0.9999*mn*sin(x1)*sin(x1);
Vw=x2*x2;
Vr=x3;
mub=x4;
mue=x5;
//proton:
m8= mp-xs*Vs; //massa efetiva
k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Neutron
m8= mn-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Lambda
m8= ml-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Sigma -
m8= msm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Sigma +
m8= msp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Sigma 0
m8= ms0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Chi -
m8= mcm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
//Chi 0
m8= mc0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
rho_s+=dens_s(m8,k);
}
f=b*mn*Vs*Vs;
f=f+c*pow(Vs,3.0);
f=f-xs*rho_s;
f=Vs+gms2*f;
return f;
}
double Eq2 (double x1, double x2, double x3, double x4, double x5, void *params) //equacao para omega
{
double gmw2 = ((struct rparams *) params)->gmw2;
double xw = ((struct rparams *) params)->xw;
double xs = ((struct rparams *) params)->xs;
double Vs, Vr, Vw, mub, mue;
double n=0.0;
double f,m8, k2,k;
//Mapeio:
Vs=0.9999*mn*sin(x1)*sin(x1);
Vw=x2*x2;
Vr=x3;
mub=x4;
mue=x5;
//proton:
m8= mp-xs*Vs; //massa efetiva
k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Neutron
m8= mn-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Lambda
m8= ml-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma -
m8= msm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma +
m8= msp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma 0
m8= ms0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi -
m8= mcm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi 0
m8= mc0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
f=(xw*gmw2)*n;
f=Vw-f;
return f;
}
double Eq3 (double x1, double x2, double x3, double x4, double x5, void *params) //equacao para o rho
{
double gmr2 = ((struct rparams *) params)->gmr2;
double xs = ((struct rparams *) params)->xs;
double xr = ((struct rparams *) params)->xr;
double Vs, Vr, Vw, mub, mue;
double n=0.0;
double f,m8,k2,k;
//Mapeio:
Vs=0.9999*mn*sin(x1)*sin(x1);
Vw=x2*x2;
Vr=x3;
mub=x4;
mue=x5;
//proton:
m8= mp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=0.5*dens(k);
}
//Neutron
m8= mn-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=-0.5*dens(k);
}
//Sigma -
m8= msm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=-1.0*dens(k);
}
//Sigma +
m8= msp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi -
m8= mcm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=-0.5*dens(k);
}
//Chi 0
m8= mc0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=0.5*dens(k);
}
f=n*xr*gmr2;
f=f-Vr;
return f;
}
double Eq4 (double x1, double x2, double x3, double x4, double x5, void *params) //neutralidade da carga eletrica (sum(ni*Qi)=0)
{
double xs = ((struct rparams *) params)->xs;
double Vs, Vr, Vw, mub, mue;
double n=0.0, ke=0.0, kmu=0.0;
double m8, k2, k;
//Mapeio:
Vs=0.9999*mn*sin(x1)*sin(x1);
Vw=x2*x2;
Vr=x3;
mub=x4;
mue=x5;
//barions
//proton:
m8= mp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma -
m8= msm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=-1.0*dens(k);
}
//Sigma +
m8= msp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi -
m8= mcm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=-1.0*dens(k);
}
//leptons
ke=mue*mue-me*me;
if(ke>0.0){
ke=sqrt(ke);
n+=-1.0*dens(ke);
}
kmu=mue*mue-mmu*mmu;
if(kmu>0.0){
kmu=sqrt(kmu);
n+=-1.0*dens(kmu);
}
return n;
}
double Eq5 (double x1, double x2, double x3, double x4, double x5,void *params) //densidade (de número) barionica (MeV^3)
{
double xs = ((struct rparams *) params)->xs;
double rho = ((struct rparams *) params)->rho;
double Vs, Vr, Vw, mub, mue;
double n=0.0;
double f,m8,k2,k;
//Mapeio:
Vs=0.9999*mn*sin(x1)*sin(x1);
Vw=x2*x2;
Vr=x3;
mub=x4;
mue=x5;
//proton:
m8= mp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Neutron
m8= mn-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Lambda
m8= ml-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma -
m8= msm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma +
m8= msp-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,1.0,1.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Sigma 0
m8= ms0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.0,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi -
m8= mcm-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,-1.0,-0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
//Chi 0
m8= mc0-xs*Vs;
k2=kf2(mub,mue,m8,Vw,Vr,0.0,0.5,params); //qe, I3
if(k2>0.0){
k=sqrt(k2);
n+=dens(k);
}
f=n-rho;
return f;
}
double kf2(double mub, double mue, double m8, double Vw, double Vr, double qe, double I3, void *params)
{
double xr = ((struct rparams *) params)->xr;
double xw = ((struct rparams *) params)->xw;
double k2, mu;
mu=mub-qe*mue;
k2=mu-xw*Vw-I3*xr*Vr; //momento de Fermi ao quadrado
k2=k2*k2;
k2=k2-m8*m8;
// printf("k2=%9.9lf\n",k2);
return k2;
}
double dens_s(double m8, double k)
{
double Pi=M_PI;
double ns;
ns=k+sqrt(k*k+m8*m8);
ns=m8*m8*log(ns/m8);
ns=k*sqrt(k*k+m8*m8)-ns;
ns=m8/(2.0*Pi*Pi)*ns;
return ns;
}
double dens(double k)
{
double Pi=M_PI;
double n;
n=pow(k,3.0);
n=n/(3.0*Pi*Pi);
return n;
}
