- #1
Gaso
- 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:
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. :
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 ,[/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.
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;
}
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";
}
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 ,[/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: