libode
Easy to compile, fast ODE integrators as C++ classes
Loading...
Searching...
No Matches
ode_lobatto_iiic_6.cc
Go to the documentation of this file.
1
2
4
5namespace ode {
6
7void NewtonLobattoIIIC6::f_Newton (double *x, double *y) {
8
9 (void)x; //supress unused variable warning
10 unsigned long i;
11 int j,n;
12 double dt = *dt_; //get time step from solver object
13
14 for (n=0; n<nk_; n++) {
15 //evaluate temporary solution based on current k values
16 for (i=0; i<neq_; i++) {
17 soltemp_[i] = sol_[i];
18 for (j=0; j<nk_; j++) soltemp_[i] += dt*a[n][j]*k_[j][i];
19 }
20 //compute k1 values
22 //evaluate the Newton function's k1 component
23 for (i=0; i<neq_; i++) y[i+n*neq_] = k_[n][i] - ftemp_[i];
24 }
25}
26
27void NewtonLobattoIIIC6::J_Newton (double *x, double **J) {
28
29 (void)x; //supress unused variable warning
30 unsigned long i,j;
31 int m,n;
32 double dt = *dt_; //get time step from solver object
33
34 if ( get_modified() ) jac(sol_, Jac_);
35
36 for (n=0; n<nk_; n++) {
37 if ( !get_modified() ) {
38 //evaluate temporary solution based on current k values for k1
39 for (i=0; i<neq_; i++) {
40 soltemp_[i] = sol_[i];
41 for (m=0; m<nk_; m++) soltemp_[i] += dt*a[n][m]*k_[m][i];
42 }
43 //evaluate J at soltemp
45 }
46 //evaluate the Newton jacobian
47 for (i=0; i<neq_; i++) {
48 for (j=0; j<neq_; j++) {
49 for (m=0; m<nk_; m++) {
50 J[i+n*neq_][j+m*neq_] = -dt*a[n][m]*Jac_[i][j];
51 }
52 }
53 J[i+n*neq_][i+n*neq_] += 1.0; //subtract from the identity matrix
54 }
55 }
56}
57
58//------------------------------------------------------------------------------
59
61 OdeAdaptive (neq, true),
62 OdeIRK (neq, 4) {
63
64 method_ = "LobattoIIIC6";
65
66 int nk = 4;
67 a = new double*[nk];
68 for (int i=0; i<nk; i++) a[i] = new double[nk];
69 b = new double[nk];
70
71 double r = sqrt(5.0);
72
73 a[0][0] = 1.0/12; a[0][1] = -r/12; a[0][2] = r/12; a[0][3] = -1.0/12;
74 a[1][0] = 1.0/12; a[1][1] = 1.0/4; a[1][2] = (10 - 7*r)/60; a[1][3] = r/60;
75 a[2][0] = 1.0/12; a[2][1] = (10 + 7*r)/60; a[2][2] = 1.0/4; a[2][3] = -r/60;
76 a[3][0] = 1.0/12; a[3][1] = 5.0/12; a[3][2] = 5.0/12; a[3][3] = 1.0/12;
77
78 b[0] = 1.0/12; b[1] = 5.0/12; b[2] = 5.0/12; b[3] = 1.0/12;
79
80 newton_ = new NewtonLobattoIIIC6(neq, nk, this);
81 newton_->set_modified(true);
82}
83
85 for (int i=0; i<nk_; i++) delete [] a[i];
86 delete [] a;
87 delete [] b;
88}
89
90void OdeLobattoIIIC6::step_ (double dt) {
91
92 unsigned long i;
93 //set k values to zero
94 for (i=0; i<neq_*nk_; i++) kall_[i] = 0.0;
95 //solve for slopes
96 newton_->solve_Newton(kall_);
97 //compute new solution
98 for (i=0; i<neq_; i++) sol_[i] += dt*(b[0]*k_[0][i]
99 + b[1]*k_[1][i]
100 + b[2]*k_[2][i]
101 + b[3]*k_[3][i]);
102}
103
104} // namespace ode
Nonlinear system solver for OdeLobattoIIIC6.
Base class implementing solver functions with adaptive time steps.
unsigned long neq_
number of equations in the system of ODEs
Definition ode_base.h:327
std::string method_
the "name" of the solver/method, as in "Euler" or "RK4"
Definition ode_base.h:319
double * sol_
array for the solution, changing over time
Definition ode_base.h:333
Provides a large vector containing the slope values of all stages with pointers to each of the indivi...
Definition ode_irk.h:9
double ** k_
individual k arrays for each stage
Definition ode_irk.h:25
int nk_
number of RK stages
Definition ode_irk.h:21
double * kall_
pointer to single array storing all stage k values
Definition ode_irk.h:23
OdeLobattoIIIC6(unsigned long neq)
constructs
double ** Jac_
pointer to the solver's Jacobian matrix
void fun(double *solin, double *fout)
wrapper around system evaluation function
unsigned long neq_
ODE system size.
double * dt_
pointer to time step member
void jac(double *solin, double **Jout)
wrapper around Jacobian evaluation function
double * ftemp_
temporary values for evaluation of Newton function
double * sol_
pointer to the solver's solution vector
double * soltemp_
temporary solution values
void set_modified(bool modified)
sets whether modified Newtion's is being used
Definition ode_newton.h:56
int solve_Newton(double *x)
Solve the system of equations.
Definition ode_newton.cc:84
bool get_modified()
gets whether modified Newtion's is being used
Definition ode_newton.h:45
int nk_
number of stages or k vectors
double ** a
pointer to tableau coefficients
double ** k_
pointer to the stage slopes of RK methods