Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Complex ODE system numerically - GSL ODE SOLVER

  1. Mar 19, 2009 #1

    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 (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,
               if (status != GSL_SUCCESS)
               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: May 4, 2017
  2. jcsd
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