# Complex ODE system numerically - GSL ODE SOLVER

1. Mar 19, 2009

### 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:
Code (Text):
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 (Text):

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 [Broken],[/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.

Last edited by a moderator: May 4, 2017
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook

Can you offer guidance or do you also need help?
Draft saved Draft deleted