Complex ODE system numerically - GSL ODE SOLVER

  • Thread starter Gaso
  • Start date
4
0
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:

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

[tex]E_{mn} = E_{m} - E_{n}[/tex], where [tex]E_{n}[/tex] are the eigenvalues.

[tex]E_{n}[/tex], [tex]\Psi_{n}[/tex] are Real, but [tex]C_{n}[/tex] 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 [tex]C_{n}[/tex] are very large, but I think that the values should be [tex]\left| C_{n} \right| \leq 1[/tex], so i suspect that my ODE system in GSL isn't well defined in func and jac.

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

[tex]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} [/tex]
[tex] 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}[/tex]; 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 [tex]C_{n}(0) = 0[/tex], except for the Eigenstate k, which is the state of the system at t=0, [tex]C_{k}(0) = 1[/tex], all imaginary parts of C at t=0 are ZERO.

Function func in my program for test case N=8:
Code:
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 Im part EigValues, vrr - real part Eigen Vectors, vri = 0, imaginary part Eigen Vectors. :
Code:
       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 [tex]C_{k}(t)[/tex] calculated with Mathematica, for my time dependet potential, and Ck is always under 1:
http://www.shrani.si/f/3k/VB/3NAQF9A9/cnexpected.jpg [Broken],[/URL]
but GSL ODE solver gives me the value of [tex]C_{k}(t) \rightarrow \infty [/tex], for example [tex]C_{k}(8.4) \approx 10^{7} [/tex]

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.
 
Last edited by a moderator:

Related Threads for: Complex ODE system numerically - GSL ODE SOLVER

Replies
0
Views
1K
Replies
4
Views
7K
Replies
5
Views
761
Replies
3
Views
2K
Replies
6
Views
9K
Replies
0
Views
2K
  • Last Post
Replies
10
Views
1K
  • Last Post
Replies
2
Views
2K

Hot Threads

Top