Runga-Kutta


Bill Casey
09-30-2003, 08:23 AM
Folks-

Does anyone have a simple program to share driving Num Recipes that can handle a system of nonlinear first-order ODEs? These are equations for a demonstration of chemical kinetics.

I am trying to save the time of debugging. This application is so common that someone must have already written one.

Bill

movax
07-03-2005, 04:56 PM
Here is onp library for numerical and symbolic computation:

http://tsk.ch.uj.edu.pl/cgi-bin/viewcvs.cgi/onp.tar.gz?view=tar

in directory /onp/src/onp/onp_math/ode/ there is ODEs related routinges.

(Or browse manualy:
http://tsk.ch.uj.edu.pl/cgi-bin/viewcvs.cgi/onp/ )

These files are under BSD licencse, You can use them for any purpose. You can easy get only part which computeate ODEs.

This methods was succesfully used to simulate nonlinear dynamics in chemical and physical systems (implementatnio and example code that use onp code here: http://tsk.ch.uj.edu.pl/cgi-bin/viewcvs.cgi/chaos/ in directory systems/ there is definitions of differential equations systems, and i.e. in ode.c there call to misc ode solving methods from onp)

The adaptative step size Kash-Carp method (in rkf5_ck.c) for system of equations is implemented by me, and can contain errors, or small differences in comparision to this in NR, (of course all coefficents are the same), but it works very good.

Witek

res0gb72
08-31-2005, 11:16 AM
Ok, so I have a series of differential equations which model a chemical reaction. I need an ODE solver, I have never used C++. I have Visual Studio 2003 and the DEV C++. I downloaded the file that you linked but I don't know what to do next. Can you please explain exact steps for beginners. Thank you.

movax
09-02-2005, 04:17 AM
Last time there was changes in onp API, but you probably downloaded old version.

You need entire ode/ directory from there.

here is example program in old API.

#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#include <ode.h>


/*
dx/dt = 0.52*x - 0.13*y
dy/dt = 5.77*x - 0.5*y
*/

/* in old API, for each equation in ode system we must have separate C function */
/* n is size of system+1, here 3, i is number of equation, in our example, 2 or 3, x0 is variables */
ode_real example_f(int n, int i, ode_real *x0) {
switch (i) {
case 1: return 0.52*x0[1] - 0.13*x0[2];
case 2: return 5.77*x0[1] - 0.5*x0[2];
}
}

/* then we need have array for variables */
ode_real* example_var[3];



int main(int argc, char **argv) {
/* initialization */
example_var[0] = 0.0; /* time */
example_var[0] = 0.1; /* x */
example_var[0] = 0.0; /* y */
/* run integrator */
/*rk4_integrator(int n, f_pars f, ode_real *x, ode_real delta_t, ode_real epsilon_y, ode_real time, ode_real delta_t_view)*/
rk4_integrator(3, example_f, example_var, 0.01, 0.01, 10.0, 0.1);
return 0;
}

/* this program will print to standard output evolution of this equation from t=0, to t=10, stepsize=0.01, printing will be done in 0.1 intervals. */

in new api (which is more powerfull, mainly for big systems


#include <stdlib.h>
#include <stdio.h>
#include <math.h>

/* we need to include all onp procedures */
#include <onp/onp.h>

/*
dx/dt = 0.52*x - 0.13*y
dy/dt = 5.77*x - 0.5*y
*/

ode_real example_f(int n, int i, ode_real *x0, ode_real *res, void *user_ptr) {
res[1] = 0.52*x0[1] - 0.13*x0[2];
res[2] = 5.77*x0[1] - 0.5*x0[2];
return 0;
}

/* then we need have array for variables */
ode_real* example_var[3];

int main(int argc, char **argv) {
/* initialization */
example_var[0] = 0.0; /* time */
example_var[0] = 0.1; /* x */
example_var[0] = 0.0; /* y */
/* run integrator */
/*int rk4_integrator(int n, f_pars2 f, ode_real *x, ode_real delta_t,
ode_real epsilon_y, ode_real time, ode_real delta_t_view,
f_rap raportf, f_pars1 checkf, void *user_ptr)
*/
rk4_integrator(3, example_f, example_var, 0.01, 0.01, 10.0, 0.1, ode_out1, NULL);
return 0;
}


/* this do the same, but have higher performance. I have situation when i have system with about 100 equations, new API was 10 times quicker in computations. */
/* you can make your own output function (raportf)
ie.
int example_out(int n, f_pars2 f, ode_real *x, f_pars1 checkf) {
.... save array x to file, or plot it in real time, or do something other...
}
*/
/* you can also provide checkf, which should be constant during real solution, for example energy of system, when integrator will detect big error between initial energy and current, it will stop.
*/


hevaing some of this file in example.c, here is way to compile it (gcc 3.x, Linux):

#first we compile method
cd ~/ONPDIR/onp/src/onp/onp_math/ode
make
#we should have now file ode.o with all procedures
#now compile our program and link it agains ode.o
cd EXAMPLEDIR
gcc -Wall -O2 -I~/ONPDIR/onp/src/onp/onp_math/ode ~/ONPDIR/onp/src -o example example.c ~/ONPDIR/onp/src/onp/onp_math/ode/ode.o -lm


#run
./example

or
./example > example.txt
gnuplot
>plot "example.txt" u 2:3 w d
# it should be kind of spiral.


i dont know how to compile this in Visual Studio, but it is possible and i suppose not very hard.

res0gb72
09-02-2005, 11:37 AM
So I copied and pasted the new API, and got this from the compile Log:
Compiler: Default compiler
Compiler: Default compiler
Executing g++.exe...
g++.exe "C:\Dev-Cpp\Untitled1.cpp" -o "C:\Dev-Cpp\Untitled1.exe" -g3 -I"C:\Dev-Cpp\lib\gcc\mingw32\3.4.2\include" -I"C:\Dev-Cpp\include\c++\3.4.2\backward" -I"C:\Dev-Cpp\include\c++\3.4.2\mingw32" -I"C:\Dev-Cpp\include\c++\3.4.2" -I"C:\Dev-Cpp\include" -L"C:\Dev-Cpp\lib" -g3
C:\Dev-Cpp\Untitled1.cpp:6:21: onp/onp.h: No such file or directory
C:\Dev-Cpp\Untitled1.cpp:13: error: `ode_real' does not name a type
C:\Dev-Cpp\Untitled1.cpp:20: error: expected constructor, destructor, or type conversion before '*' token
C:\Dev-Cpp\Untitled1.cpp:20: error: expected `,' or `;' before '*' token
C:\Dev-Cpp\Untitled1.cpp: In function `int main(int, char**)':
C:\Dev-Cpp\Untitled1.cpp:24: error: `example_var' undeclared (first use this function)
C:\Dev-Cpp\Untitled1.cpp:24: error: (Each undeclared identifier is reported only once for each function it appears in.)

C:\Dev-Cpp\Untitled1.cpp:32: error: `example_f' undeclared (first use this function)
C:\Dev-Cpp\Untitled1.cpp:32: error: `ode_out1' undeclared (first use this function)
C:\Dev-Cpp\Untitled1.cpp:32: error: `rk4_integrator' undeclared (first use this function)

Execution terminated


I extracted the library into the Dev-C++ folder an dit still can't find it.

movax
09-03-2005, 01:24 AM
Sorry, new API is not published yet, im testing it now. I will upload all files in next week. :)

Use example for old. Error you have is something with your paths to headers file; include <ode.h> not <onp/onp.h>.
Or "ode.h" and have all files in the same dir that your files are.

res0gb72
09-06-2005, 09:42 AM
error:
Compiler: Default compiler
Compiler: Default compiler
Executing g++.exe...
g++.exe "C:\Dev-Cpp\test.cpp" -o "C:\Dev-Cpp\test.exe" -I"C:\Dev-Cpp\lib\gcc\mingw32\3.4.2\include" -I"C:\Dev-Cpp\include\c++\3.4.2\backward" -I"C:\Dev-Cpp\include\c++\3.4.2\mingw32" -I"C:\Dev-Cpp\include\c++\3.4.2" -I"C:\Dev-Cpp\include" -L"C:\Dev-Cpp\lib"
C:\Dev-Cpp\test.cpp: In function `int main(int, char**)':
C:\Dev-Cpp\test.cpp:29: error: cannot convert `double' to `ode_real*' in assignment

C:\Dev-Cpp\test.cpp:30: error: cannot convert `double' to `ode_real*' in assignment
C:\Dev-Cpp\test.cpp:31: error: cannot convert `double' to `ode_real*' in assignment

C:\Dev-Cpp\test.cpp:34: error: cannot convert `ode_real**' to `ode_real*' for argument `3' to `int rk4_integrator(int, ode_real (*)(int, int, ode_real*), ode_real*, ode_real, ode_real, ode_real, ode_real)'

Execution terminated

movax
09-08-2005, 07:56 AM
Should be:
ode_real example_var[3];

without asterix after ode_real.

Also there is error in initialization. Should be:
/* initialization */
example_var[0] = 0.0; /* time */
example_var[1] = 0.1; /* x */
example_var[2] = 0.0; /* y */


I was writing this example directly to forum :) so i didn't test it.

res0gb72
09-08-2005, 04:24 PM
Compiler: Default compiler
Executing g++.exe...
g++.exe "C:\Dev-Cpp\test.cpp" -o "C:\Dev-Cpp\test.exe" -I"C:\Dev-Cpp\lib\gcc\mingw32\3.4.2\include" -I"C:\Dev-Cpp\include\c++\3.4.2\backward" -I"C:\Dev-Cpp\include\c++\3.4.2\mingw32" -I"C:\Dev-Cpp\include\c++\3.4.2" -I"C:\Dev-Cpp\include" -L"C:\Dev-Cpp\lib"
C:\DOCUME~1\Maksim\LOCALS~1\Temp/ccoPbaaa.o(.text+0xfc):test.cpp: undefined reference to `rk4_integrator(int, double (*)(int, int, double*), double*, double, double, double, double)'
collect2: ld returned 1 exit status

Execution terminated

movax
09-09-2005, 01:47 AM
OK, as i see compilation of example.c is succesfull, but you must now link it agaist ode.o

i don know how to do this in dev-cpp, so i give you some brief manual

cd C:\Dev-Cpp
g++ -c -o test.o test.cpp -I"C:\Dev-Cpp\lib\gcc\mingw32\3.4.2\include" -I"C:\Dev-Cpp\include\c++\3.4.2\backward" -I"C:\Dev-Cpp\include\c++\3.4.2\mingw32" -I"C:\Dev-Cpp\include\c++\3.4.2" -I"C:\Dev-Cpp\include" -L"C:\Dev-Cpp\lib"

cd C:\Dev-Cpp\ode\
cd methods
gcc -c -o rk4.o rk4.c -I..\
gcc -c -o euler.o euler.c -I..\
...
gcc -c -o rkf5_ck.o rkf5_ck.c -I..\
ld -r -o methods.o rk4.o ... rkf5_ck.o
cd ..
gcc -c -o ode.o ode.c
...
ld -r -o onp_ode.o ode.c ... methods\methods.o


now link onp_ode.o agaist test.o (and never name program "test", on unix there is such command)

ld -o example.exe test.o ode\onp_ode.o

i you have any other problems you must learn how gcc and dev-cpp works, and i will not help you.
On unix you only have to execute "make" command, in Makefile there are command that will compile everything and link.