Re: Help! Questions for the developers of GSL ODE: Complex-valued ODEs and the "jac" function in GSL ODE solver?

[eknecronzontas <eknecronzontas@yahoo.com>]

Yes GSL's ode solver only handles real equations. However, there is no
factor of two slowdown, as you're only explicitly doing with GSL what
matlab must do internally. Whenever you write two complex-valued ODE's
(no matter what software you use) you are doing nothing more than
implicitly writing 4 real ODE's. The "jac" function is the Jacobian.
For yprime_i=f_i(y_j) the Jacbian gives the matrix which is (d f_i)/(dy_j). 
I don't quite remember the ordering (whether i or j is rows or columns)
though.

The specification of the Jacobian is indeed documented at:
http://www.gnu.org/software/gsl/manual/html_node/Defining-the-ODE-System.html

I assume your boundary conditions couple these ODE's which
could be solved separately since 'x' does not appear in the
first equation.

Take care,
Andrew

----- Original Message ----
From: Michael <comtech.usa@gmail.com>
To: gsl-discuss@sources.redhat.com
Sent: Sunday, June 3, 2007 4:18:17 AM
Subject: Help! Questions for the developers of GSL ODE: Complex-valued ODEs and the "jac" function in GSL ODE solver?

Help! Complex-valued ODEs and the "jac" function in GSL ODE solver?

Hi all,

I am new to GSL ode, I started looking it because I really need to
speed my current Matlab program, which runs ODE solvers repeatedly for
thousands to millions of times. I found the GSL example ode program
difficult to understand, and here are my questions:

(1)

I have to solve the following complex-valued non-linear ODE numerically,
using GSL's ode solver. But it seems that GSL's ode solver only supports
real-valued ODE...

My ODEs are:

y' = c1 * y + c2 + c3*exp(c4* y + c5*i)
x' = c6 * y

here c1, c2, c3, c4, c5 and c6 are real numbers , and "i" is the unit of
imaginary numbers.

I am planning to decompose y into yr and yi, x into xr and xi, the real
parts and the imaginary parts. And solve separately:


yr' = c1 * yr + c2 + c3*exp(c4* yr)*cos(c4* yi + c5)
yi' = c1 * yi + c2 + c3*exp(c4* yr)*sin(c4* yi + c5)
xr' = c6 * yr
xi' = c6 * yi

--------------------

Am I right? Is this the best approach to handle the ODEs with a solver only
supports real values?

My original solutions were already too slow, now by having to solve 4
equations, it is even slower, and double the computing time...

Any better approaches?


(2)

In the sample GSL ode code, there is a function called "jac", how do I
define the "jac" function for my ODE equations, as shown above?

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  double mu = *(double *)params;
  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, 2, 2);
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

---------------------

There is no documentation talking about this "jac" function. I really
don't know how to define my "jac" function for my own ode equations.

Could you please help me?

Thanks a lot!

Mike