help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [Help-gsl] Help! Complex-valued ODEs and the "jac" function in GSL O


From: Warren Weckesser
Subject: Re: [Help-gsl] Help! Complex-valued ODEs and the "jac" function in GSL ODE solver?
Date: Sun, 03 Jun 2007 18:22:09 -0400

Hi Mike,

On Sun, 2007-06-03 at 01:09 -0700, Michael wrote:
> 
> 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...

Decomposing the complex numbers into real and imaginary parts seems
like the most straight-forward approach.  I doubt that the computing
time will be doubled, since adding and multiplying real numbers is
faster than adding and multplying complex numbers.

By the way, if c2 is real, it should not show up in the
equation for yi'.

> 
> 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?
<snip>
> 
> 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.

"jac" refers to the Jacobian matrix of the right-hand side of your
differential equations. See, for example, 
    http://en.wikipedia.org/wiki/Jacobian
or
    http://mathworld.wolfram.com/Jacobian.html

> Could you please help me?

You might be interested in trying VFGEN, the program that I mentioned
in my email to the mailing list a little while ago.  The representation
of your system in VFGEN's input format is:

=== system.vf ==========
<?xml version="1.0" ?>
<VectorField Name="system">
<Parameter Name="c1" DefaultValue="1.0" />
<Parameter Name="c2" DefaultValue="1.0" />
<Parameter Name="c3" DefaultValue="1.0" />
<Parameter Name="c4" DefaultValue="1.0" />
<Parameter Name="c5" DefaultValue="1.0" />
<Parameter Name="c6" DefaultValue="1.0" />
<StateVariable
    Name="yr"
    Formula="c1*yr + c2 + c3*exp(c4*yr)*cos(c4*yi + c5)" />
<StateVariable
    Name="yi"
    Formula="c1*yi + c3*exp(c4*yr)*sin(c4*yi + c5)" />
<StateVariable
    Name="xr"
    Formula="c6*yr" />
<StateVariable
    Name="xi"
    Formula="c6*yi" />
</VectorField>
=== end of system.vf ====

The command

$ vfgen gsl:demo=yes system.vf

produces four files: system_gvf.c, system_gvf.h, system_solve.c, and
Makefile-system.  The C code in system_gvf.c defines the vector field
(i.e. the right-hand side) and the Jacobian. (It also defines the
derivatives of the vector field with respect to the parameters, but
you can ignore that function.)  Prototypes for the functions are in
system_gvf.h.  The file system_solve.c is a "demonstration" solver for
your system.  Compile it with the command

$ make -f Makefile-system

and run it like this:

$ ./system_solve yr=1.0 yi=0.0 xr=0.0 xi=-1.0 abserr=1e-10 > out.dat

This will generate a solution to the system; out.dat will contain
the solution, with one line for each point computed.  The format is
   t  yr  yi  xr  xi
You can also specify the parameter values, e.g.

$ ./system_solve yr=1.0 xr=-1.0 c1=3.0 c4=0.04 abserr=1e-10 > out.dat

This solver might be useful as it is, but another reason for creating
the demonstration program is to see an example of how to call the
GSL library functions to solve the system.

You can get VFGEN here:
    http://math.colgate.edu/~wweckesser/software/vfgen


Best regards,

Warren Weckesser





reply via email to

[Prev in Thread] Current Thread [Next in Thread]