[skabilan]

Hi,

I am new to OpenFoam and was trying to use the ODE
solver. I took a look at the ODETest but was not able to understand how the equations are cast. can some one show me how I can solve a simple system like the following?

C --- INITIAL VALUES
X=0.0D0
Y(1)=2.0D0
Y(2)=-0.66D0
C --- ENDPOINT OF INTEGRATION
XEND=2.0D0
XDELTA=0.2
C --- SYSTEM OF EQUATIONS
F(1)=Y(2)
F(2)=((1-Y(1)**2)*Y(2)-Y(1))/2.0


ODETest
\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "IOmanip.H"
#include "ODE.H"
#include "ODESolver.H"
#include "RK.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

class testODE
:
public ODE
{

public:

testODE()
{}

label nEqns() const
{
return 4;
}

void derivatives
(
const scalar x,
const scalarField& y,
scalarField& dydx
) const
{
dydx[0] = -y[1];
dydx[1] = y[0] - (1.0/x)*y[1];
dydx[2] = y[1] - (2.0/x)*y[2];
dydx[3] = y[2] - (3.0/x)*y[3];
}

void jacobian
(
const scalar x,
const scalarField& y,
scalarField& dfdx,
Matrix<scalar>& dfdy
) const
{
dfdx[0] = 0.0;
dfdx[1] = (1.0/sqr(x))*y[1];
dfdx[2] = (2.0/sqr(x))*y[2];
dfdx[3] = (3.0/sqr(x))*y[3];

dfdy[0][0] = 0.0;
dfdy[0][1] = -1.0;
dfdy[0][2] = 0.0;
dfdy[0][3] = 0.0;

dfdy[1][0] = 1.0;
dfdy[1][1] = -1.0/x;
dfdy[1][2] = 0.0;
dfdy[1][3] = 0.0;

dfdy[2][0] = 0.0;
dfdy[2][1] = 1.0;
dfdy[2][2] = -2.0/x;
dfdy[2][3] = 0.0;

dfdy[3][0] = 0.0;
dfdy[3][1] = 0.0;
dfdy[3][2] = 1.0;
dfdy[3][3] = -3.0/x;
}
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// Main program:

int main(int argc, char *argv[])
{
argList::validArgs.clear();
argList::validArgs.append("ODESolver");
argList args(argc, argv);

word ODESolverName(args.args()[1]);

testODE ode;
autoPtr<odesolver> odeSolver = ODESolver::New(ODESolverName, ode);

scalar xStart = 1.0;
scalarField yStart(ode.nEqns());
yStart[0] = ::Foam::j0(xStart);
yStart[1] = ::Foam::j1(xStart);
yStart[2] = ::Foam::jn(2, xStart);
yStart[3] = ::Foam::jn(3, xStart);

scalarField dyStart(ode.nEqns());
ode.derivatives(xStart, yStart, dyStart);

Info<< setw(10) << "eps" << setw(12) << "hEst";
Info<< setw(13) << "hDid" << setw(14) << "hNext" << endl;
Info<< setprecision(6);

for (label i=0; i<15; i++)
{
scalar eps = ::Foam::exp(-scalar(i + 1));

scalar x = xStart;
scalarField y = yStart;
scalarField dydx = dyStart;

scalarField yScale(ode.nEqns(), 1.0);
scalar hEst = 0.6;
scalar hDid, hNext;
odeSolver->solve(ode, x, y, dydx, eps, yScale, hEst, hDid, hNext);

Info<< scientific << setw(13) << eps;
Info<< fixed << setw(11) << hEst;
Info<< setw(13) << hDid << setw(13) << hNext
<< setw(13) << y[0] << setw(13) << y[1]
<< setw(13) << y[2] << setw(13) << y[3]
<< endl;
}

scalar x = xStart;
scalar xEnd = x + 1.0;
scalarField y = yStart;

scalarField yEnd(ode.nEqns());
yEnd[0] = ::Foam::j0(xEnd);
yEnd[1] = ::Foam::j1(xEnd);
yEnd[2] = ::Foam::jn(2, xEnd);
yEnd[3] = ::Foam::jn(3, xEnd);

scalar hEst = 0.5;

odeSolver->solve(ode, x, xEnd, y, 1e-4, hEst);
Info<< nl << "Analytical: y(2.0) = " << yEnd << endl;
Info << "Numerical: y(2.0) = " << y << ", hEst = " << hEst << endl;

Info << "\nEnd\n" << endl;

return 0;
}


// ************************************************** *********************** //