Resolution of Peeble's ODE equation

  • Context: Undergrad 
  • Thread starter Thread starter redgun
  • Start date Start date
  • Tags Tags
    Recombination Universe
Click For Summary
SUMMARY

The forum discussion centers on resolving Peebles' ODE equation for modeling the recombination history of the universe using C programming. The user encounters divergence in the electron fraction, Xe, when approaching the Saha limit due to excessively large values on the right-hand side of the ODE. The code provided utilizes physical constants and numerical methods, including the Runge-Kutta method, to compute electron density and fractions. The user seeks assistance in identifying potential issues within their implementation that prevent convergence of the results.

PREREQUISITES
  • Understanding of ODE (Ordinary Differential Equations) and numerical methods, specifically Runge-Kutta.
  • Familiarity with cosmological parameters and their physical significance, such as OmegaB0 and H_0.
  • Proficiency in C programming, particularly in memory management and function pointers.
  • Knowledge of physical constants relevant to cosmology, such as the speed of light (c) and electron mass (m_e).
NEXT STEPS
  • Investigate GSL (GNU Scientific Library) methods for improved stability in solving ODEs.
  • Explore the implications of the Saha limit on electron fraction calculations in cosmological models.
  • Learn about spline interpolation techniques for better handling of numerical data in C.
  • Examine the impact of varying initial conditions on the convergence of ODE solutions.
USEFUL FOR

Researchers in cosmology, software developers working on astrophysical simulations, and anyone interested in numerical methods for solving differential equations in a scientific context.

redgun
Messages
2
Reaction score
0
TL;DR
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){
count << "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;
count << "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));

count << "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];
count << q << endl;
double a=exp(x);
double n_b = (OmegaB0*3*H_0*H_0)/(8*M_PI*G*m_h*a*a*a);
count << n_b << endl;
double n1s = (1-q)*n_b;
count << 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;
count << Xe << endl;
ne_array = ne;
}
for (int i = 0; i < npts + 1; i++) {
count << "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.
 
Last edited:
Space news on Phys.org
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);

count << 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 :
1742315968239.png

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 ?
 

Similar threads

C/C++ Why Am I Getting Linker Errors When Using Const Methods in C++?
  • · Replies 13 ·
Replies
13
Views
2K
Graduate Solve RK4 in C++ for Sekhar Mass White Dwarf Star
  • · Replies 2 ·
Replies
2
Views
7K