I Resolution of Peeble's ODE equation

[redgun]

TL;DR Summary

Here is my code for recombination history of universe that have a problem when solving peebles ode

Hello everyone,

I am trying to code the recombination history of the universe in C. More precisely I have difficulties to solve the peebles ode equation. In fact when I am close to the saha limit, I begin the resolution but the rhs of the ode seems to be too big and it makes my result Xe diverge. Someone have an idea to solve the problem ?

Here is my code :

[CODE lang="c" title="Recombination"]#include "recombinaison.h"
using namespace std;

// Constantes physiques

const double k_b = 1.38064852e-23;
const double c = 299792458;
double const G = 6.67430e-11;
const double m_e = 9.10938356e-31;
const double e = 1.60217653e-19;
const double m_h = 1.6735575e-27;
const double h_p = (6.62607015e-34)/(2*M_PI);
const double h = 6.62607015e-34;
const double eps_0 = 13.605693122994; // en eV
const double alpha = 1/137.0359991;
const double mu_0 = 4*M_PI*1e-7;
const double e0 = 1/(mu_0*c*c);
const double khi_0 = 24.587387; // en eV
const double khi_1 = 4.0*eps_0; // en eV
const double Ry = (m_e*pow(e,4))/(8*e0*e0*h*h*e); // en eV;
const double L2s1s = 8.227;
const double sigmaT = (8*M_PI*alpha*alpha*h_p*h_p)/(3*c*c*m_e*m_e);

// Créer et initialiser Recombinaison

Recombinaison* creer_recombinaison(Cosmologie* cosmo, double Yp){
Recombinaison* rec = (Recombinaison*)malloc(sizeof(Recombinaison));
if (rec==NULL){
cout << "Erreur : allocation mémoire échouée pour Recombinaison" << endl;
return NULL;
}

rec -> cosmo = cosmo;
rec -> Yp = Yp;

rec -> x_start = -15.0;
rec -> x_end = 0.0;
rec -> npts = 10000;
rec -> Xe_saha_limit = 0.9900;

rec -> Xe_spline = NULL;
rec -> log_ne_spline = NULL;
rec -> tau_spline = NULL;
rec -> g_tilde_spline=NULL;

return rec;
}

// Fonction pour libérer l'espace

void free_recombinaison(Recombinaison* rec){
if (rec == NULL) return;
if (rec->Xe_spline != NULL){
free_spline(rec->Xe_spline);
rec->Xe_spline = NULL;
}
if (rec->tau_spline != NULL){
free_spline(rec->tau_spline);
rec->tau_spline = NULL;
}
if (rec->g_tilde_spline != NULL){
free_spline(rec->g_tilde_spline);
rec->g_tilde_spline = NULL;
}
free(rec);
}

// Calcul des fractions d'électrons


void rk4_rec(void(*f)(Recombinaison*,double*,double,double*), Recombinaison* rec, double* q, double x, double dx){
double *p1 = (double*)malloc(sizeof(double));
double *p2 = (double*)malloc(sizeof(double));
double *p3 = (double*)malloc(sizeof(double));
double *p4 = (double*)malloc(sizeof(double));
double *q_int = (double*)malloc(sizeof(double));

f(rec,q,x,p1);
*q_int = *q + (dx/2.)*(*p1);
f(rec,q_int,x+dx/2.,p2);
*q_int = *q + (dx/2.)*(*p2);
f(rec,q_int,x+dx/2.,p3);
*q_int = *q + dx*(*p3);
f(rec,q_int,x+dx,p4);
*q = *q + (dx/6.)*(*p1+2.*(*p2)+2.*(*p3)+*p4);
free(p1); free(p2); free(p3); free(p4); free(q_int);
}

void rk4_rec_p(void(*f)(Recombinaison*,double*,double,double*,double), Recombinaison* rec, double* q, double x, double dx, double n1s){
double *p1 = (double*)malloc(sizeof(double));
double *p2 = (double*)malloc(sizeof(double));
double *p3 = (double*)malloc(sizeof(double));
double *p4 = (double*)malloc(sizeof(double));
double *q_int = (double*)malloc(sizeof(double));

f(rec,q,x,p1,n1s);
*q_int = *q + (dx/2.)*(*p1);
f(rec,q_int,x+dx/2.,p2,n1s);
*q_int = *q + (dx/2.)*(*p2);
f(rec,q_int,x+dx/2.,p3,n1s);
*q_int = *q + dx*(*p3);
f(rec,q_int,x+dx,p4,n1s);
*q = *q + (dx/6.)*(*p1+2.*(*p2)+2.*(*p3)+*p4);
free(p1); free(p2); free(p3); free(p4); free(q_int);
}

void fraction_electron_saha(Recombinaison* rec, double x, double* Xe, double* ne) {
double a = exp(x);
double z = (1/a)-1;
const double OmegaB0 = (rec->cosmo)->OmegaB_0;
const double H_0 = get_H_0(rec->cosmo);
const double TCMB_0 = (rec->cosmo)->TCMB_0;
const double Yp = rec->Yp;
double T = TCMB_0/a;

double n_b = (3*OmegaB0*H_0*H_0)/(8*M_PI*G*m_h*a*a*a);
double h_c = get_h(rec->cosmo);

double A = 26.33+1.5*log(TCMB_0)-log(OmegaB0*h_c*h_c*pow(1+z,1.5))-157681.1594/(TCMB_0*(1+z));
double B = exp(A);

*Xe = (-B+sqrt(B*B+4.0*B))*0.5;
*ne = (*Xe)*n_b;
cout << "Xe : " << *Xe << ", T : "<< T << ", z : " << z << endl;
}


void sys_diff_peebles(Recombinaison* rec, double* q, double x, double* qp, double n1s){
double a=exp(x);
const double OmegaB0 = (rec->cosmo)->OmegaB_0;
const double H_0 = get_H_0(rec->cosmo);
const double TCMB_0 = (rec->cosmo)->TCMB_0;
double T = TCMB_0/a;
double b = (eps_0*e)/(k_b*T);
double H = H_x(rec->cosmo,x);
double n_b = (3*OmegaB0*H_0*H_0)/(8*M_PI*G*m_h*a*a*a);
double L = 8.23;


const double alpha_e = 2.6e-19*pow(T*1e-4,-0.8);
const double beta_e = (alpha_e*pow(2*M_PI*m_e*k_b*T,3./2.)*exp(-b/4.))/(h*h*h);
const double lambda_a = 121.567e-9;
const double E_a = (h*c)/(lambda_a*1.60217653e-19);
const double B = (E_a*e)/(k_b*T);
const double K = pow(lambda_a,3)/(8*M_PI*H);
double C = (1.0 + K * L* n1s) / (1.0 + K * n1s * (L+beta_e));

cout << "C : " << C << " beta_e : " << beta_e << " b : " << b << " H : " << H << " K : " << K << " n1s : " << n1s << " B : " << B << " a_e : " << alpha_e << endl;

*qp = (C/H)*(beta_e*(1-(*q))*exp(-B)-n_b*alpha_e*(*q)*(*q));
}



void solve_electron_density(Recombinaison* rec){
const int npts = rec->npts;
const double OmegaB0 = (rec->cosmo)->OmegaB_0;
const double H_0 = get_H_0(rec->cosmo);
double* x_array = (double*)malloc((npts+1)*sizeof(double));
double* Xe_array = (double*)malloc((npts+1)*sizeof(double));
double* ne_array = (double*)malloc((npts+1)*sizeof(double));

linspace(rec->x_start, rec->x_end, npts+1, x_array);
bool flag = true;
double Xe, ne;
for (int i=0;i<npts+1;i++){
double x = x_array;
if (flag){
double Xe_temp, ne_temp;
fraction_electron_saha(rec,x,&Xe_temp,&ne_temp);
Xe=Xe_temp;
ne=ne_temp;
if (Xe < rec->Xe_saha_limit) flag = false;
}
else{
double q=Xe_array[i-1];
cout << q << endl;
double a=exp(x);
double n_b = (OmegaB0*3*H_0*H_0)/(8*M_PI*G*m_h*a*a*a);
cout << n_b << endl;
double n1s = (1-q)*n_b;
cout << n1s << endl;
rk4_rec_p(sys_diff_peebles,rec,&q,x_array[i-1],x-x_array[i-1],n1s);
Xe=q;
ne=Xe*n_b;
}
Xe_array = Xe;
cout << Xe << endl;
ne_array = ne;
}
for (int i = 0; i < npts + 1; i++) {
cout << "x: " << x_array << ", Xe: " << Xe_array << endl;
}
rec->Xe_spline = calcul_spline_cubique_naturelle(x_array, Xe_array, npts+1);
rec->log_ne_spline = calcul_spline_cubique_naturelle(x_array, ne_array, npts+1);
free(x_array);
free(Xe_array);
free(ne_array);
}[/CODE]

If needed I can bring you my code for the cosmological parameters like H(x) but there are well calculated. Ask me if needed.

 

[redgun]

I decided to use gsl methods to because rk4 might cause some problems of stability. Now I try to solve Peebles ODE with this code :
[CODE lang="c" title="Peebles ode"]int peebles_ode(double x, const double y[], double f[], void* params) {
Recombinaison* rec = (Recombinaison*)params;
double a = exp(x);
const double OmegaB0 = rec->cosmo->OmegaB_0;
const double H_0 = get_H_0(rec->cosmo);
const double TCMB_0 = (rec->cosmo)->TCMB_0;
double T = TCMB_0/a;
double b = (eps_0*e)/(k_b*T);
double H = H_x(rec->cosmo, x);
double n_b = (3*OmegaB0*H_0*H_0)/(8*M_PI*G*m_h*a*a*a);

const double phi2=0.448*log(b);
const double alpha2 = (64*M_PI*alpha*alpha*h_p*h_p*sqrt(b)*phi2)/(sqrt(27*M_PI)*m_e*m_e*pow(c,4));
const double beta = alpha2*pow((m_e*k_b*T)/(2*M_PI*h_p*h_p),1.5)*exp(-b);
const double beta2 = alpha2*pow((m_e*k_b*T)/(2*M_PI*h_p*h_p),1.5)*exp(-b/4.);
const double L2s1s = 8.227;

double Xe = y[0]; // y[0] = Xe
double n1s = (1 - Xe) * n_b;
double La = (H*pow(3*eps_0*e,3))/(8*8*M_PI*M_PI*pow(h_p,3)*n1s);
double C = (L2s1s+La) / (L2s1s+La+beta2);

f[0] = (C / H) * (beta * (1 - Xe) - n_b * alpha2 * Xe * Xe);

cout << f[0] << endl;

return GSL_SUCCESS;
}[/CODE]

But there is a problem somewhere, because when I compute, it yield a result that not converges to 0 :

I have the values of f[0] that varies between 1e-13 at start and 1e-25 at the end. There is probably an issue with my ode, can someone see what is the problem here ?