Complex ODE system numerically - GSL ODE SOLVER

[Gaso]

Hi,

I need to solve very large complex ODE system. It is about time evaluation of system, which is at time t=0 in eigenstate with the smallest eigenvalue. For my test case I am trying to solve smaller similar problem, the ODE system is like:

$$C^{'}_{m} = - i \sum^{N-1}_{n=0} C_{n}(t) Exp[ i E_{mn} t] V_{mn}(t),$$
$$C_{m}(0) = 0,$$ $$m\neq k$$
$$C_{k}(0) = 1,$$
where k is index number of smallest eigenvalue(initial state)
$$V_{mn} = \langle \Psi_{m} \left| V(t) \right| \Psi_{n}\rangle$$ is a matrix element $$\Psi_{n}$$ are eigenvectors.

$$E_{mn} = E_{m} - E_{n}$$, where $$E_{n}$$ are the eigenvalues.

$$E_{n}$$, $$\Psi_{n}$$ are Real, but $$C_{n}$$ are Complex.

I tried to solve the system of size N=8 in Mathematica with NDSOLVE, and it gives me the expected result. I need to use some ODE solver in c++ to solve the system numerically, so I tried with GNU GSL ODE solver. I tried with different methods rk2, rk4, etc..., but the solutions of $$C_{n}$$ are very large, but I think that the values should be $$\left| C_{n} \right| \leq 1$$, so i suspect that my ODE system in GSL isn't well defined in func and jac.

Because the $$C_{n}$$ are Complex, I have written my system in two parts, Re and Im.

$$C^{'}_{m}(t) = \sum^{N-1}_{n=0} \left( C_{n}(t) sin(E_{mn} t) + C_{n+N}(t) cos(E_{mn} t) \right) V_{mn}$$
$$C^{'}_{m+N}(t) = \sum^{N-1}_{n=0} \left( - C_{n}(t) cos(E_{mn} t) + C_{n+N}(t) sin(E_{mn} t) \rihgt) V_{mn}$$; m=0, ..., N-1
where m= 0, ... N-1 for Real part of C, and m=N, ..., 2N-1 for Imaginary part of C.
The current system has dimension of 2N.

Because, my Eigenvalues and Eigenvectors are Real, I defined my $$C_{n}(0) = 0$$, except for the Eigenstate k, which is the state of the system at t=0, $$C_{k}(0) = 1$$, all imaginary parts of C at t=0 are ZERO.

Function func in my program for test case N=8:

int  func (double t, const double y[], double f[],
           void *params)
     {	
	SysArg *args = new SysArg[1];
 	args = *(SysArg **)params;
	vector<double>  wr = args[0].getWr();
	vector< dcovector > vrr = args[0].getVrr();
	double VMN = 0;
	double Emn = 0;
	for(long m=0; m<8; m++){
		for(long n=0; n<8; n++){
			VMN = Vmn(Potencial(t),  vrr, m, n);
			Emn = (wr[m] - wr[n]);
			f[m] +=  y[n]*sin(Emn*t)*VMN + y[n+8]*cos(Emn*t)*VMN ; 
			f[m+8] += (-y[n]*cos(Emn*t)*VMN + y[n+8]*sin(Emn*t)*VMN) ;
		}
	}
       return GSL_SUCCESS;
     }

My main program, eigenvalues and eigenvectors are calculated correctly with CPPLAPACK rutines and wr - Real part Eigenvalues, wi = 0 I am part EigValues, vrr - real part Eigen Vectors, vri = 0, imaginary part Eigen Vectors. :

       const gsl_odeiv_step_type * T 
         = gsl_odeiv_step_rk4;
     
       gsl_odeiv_step * s 
         = gsl_odeiv_step_alloc (T, 16);
         gsl_odeiv_system sys = {func, jac, 16, &arguments};
      
	 double t = 0.0, t1 = 100;
       double h = 1e-1;
       double y[16], y_err[16];
	double dydt_in[16], dydt_out[16];
// find smallest eigenvalue
	long EigVal = 0;
	double temp = wr[0];
	for(long i=1; i<vrr[0].l; i++){
    		if(wr[i] < temp){temp = wr[i]; EigVal=i;}
  	}
	
	for(int i=0; i<16; i++){
		y[i] = 0;
	}
	// EigVal - index of smallest eigenvalue, initial state
	y[EigVal] = 1;

	  /* initialise dydt_in from system parameters */
       GSL_ODEIV_FN_EVAL(&sys, t, y, dydt_in);

     while (t < t1)
         {
           int status = gsl_odeiv_step_apply (s, t, h, 
                                              y, y_err, 
                                              dydt_in, 
                                              dydt_out, 
                                              &sys);
     
           if (status != GSL_SUCCESS)
               break;
	       for(int i=0; i<16; i++){
		    dydt_in[i] = dydt_out[i];
     		}
           t += h;
		
	   myfile << t << " " << y[EigVal];
	   myfile << "\n";
         }

Runge-kuta 4 method doesn't need Jacobian, so I didn't post it here.
The correct and expected solution of time expansion of $$C_{k}(t)$$ calculated with Mathematica, for my time dependet potential, and Ck is always under 1:
http://www.shrani.si/f/3k/VB/3NAQF9A9/cnexpected.jpg ,[/URL]
but GSL ODE solver gives me the value of $$C_{k}(t) \rightarrow \infty$$, for example $$C_{k}(8.4) \approx 10^{7}$$

Do I have right approach for solving that system. Can someone check if I have defined right my system in GSL's func, there is very poor documentation for GSL ODE, and with examples included I didn't figure out, how to write the system with time dependent arguments in GSL's func and jac correctly.

 

[jvicens]

Hi,

Thank you for sharing your problem with us. Solving large and complex ODE systems can be challenging, but it seems like you have a good understanding of the problem and have approached it in a systematic way. From your code and explanation, it seems like you have defined the system correctly in GSL's func and jac. However, there are a few things that you can check to make sure that your results are correct.

Firstly, make sure that your initial conditions are correct. In your code, you have set all the imaginary parts of C to zero, except for the initial state k. Double check that this is indeed the case. Also, make sure that the initial state k is the correct one, as this can affect the results significantly.

Secondly, you can try using a different ODE solver in GSL, such as the Bulirsch-Stoer method, to see if it gives you better results. It is also worth experimenting with different step sizes and tolerances to see if that affects the results.

Lastly, it is always a good idea to compare your results with those obtained from other methods, such as Mathematica or other numerical solvers. This can help you identify any errors in your code or approach.

I hope this helps. Best of luck with solving your ODE system!