Module: Ode::Methods

Defined in:
lib/ode.rb,
lib/ode/methods.rb,
ext/ode/ode.c

Class Method Summary collapse

Class Method Details

.default_opts(name) ⇒ Object 



3
4
5
6
7
8
9
10
11
12
13
 
# File 'lib/ode/methods.rb', line 3

def self.default_opts(name)
  {
    :lsoda => {
      itol: 1, # both atol and rtol are scalar
      rtol: 1e-4,
      atol: 1e-6,
      ml: nil,
      mu: nil
    }
  }[name]
end

.lsoda(t_out, func, jac, init_t, init_y, fargs, opts) ⇒ Object 



35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
 
# File 'ext/ode/ode.c', line 35

VALUE excute_lsoda(VALUE self, VALUE t_out, VALUE func, VALUE jac, VALUE init_t, VALUE init_y, VALUE fargs, VALUE opts){
  int itol; double rtol, atol;
  int *iwork; double *rwork;
  int lrw, liw;
  int neq = RARRAY_LEN(init_y);

  double *y = (double *)ALLOC_N(double, neq);
  double t = NUM2DBL(init_t);
  double tout = NUM2DBL(t_out);

  int itask = 1;  // run solver untill t == tout
  int istate = 1; // this is the first call of the problem
  int iopt = 0;   // no optional inputs

  int jt = 2;     // jacobian type indicator. 1: user provides full jacobian. 2: interanally generated jacobian.
  VALUE ret_arr = rb_ary_new2(neq);
  double h0 = 0, hmax =0, hmin = 0;
  int ixpr=0, max_step=0, max_hnil=0, max_ordn=12, max_ords=5;
  int i;

  RESTORE_NEQ(neq)
  RESTORE_FUNC(func)
  RESTORE_FARGS(fargs)

  for(i=0; i<neq; i++){
    y[i] = NUM2DBL(rb_ary_entry(init_y, i));
  }

  // parse options (TODO: accept array)
  itol = NUM2INT(rb_hash_lookup(opts, ID2SYM(rb_intern("itol"))));

  rtol = NUM2DBL(rb_hash_aref(opts, ID2SYM(rb_intern("rtol"))));
  atol = NUM2DBL(rb_hash_aref(opts, ID2SYM(rb_intern("atol"))));

  // decide lrw, liw
  liw = 20 + neq;
  iwork = (int *)ALLOC_N(int, liw);

  if(jac != Qnil){
    jt = 1;
    RESTORE_JAC(jac);
    lrw = max(20+16*neq, 22+9*neq+neq*neq);
  }
  else{
    int ml=0, mu=0; VALUE val;
    if(val = rb_hash_aref(opts, ID2SYM(rb_intern("ml"))) != Qnil)ml = NUM2INT(val);
    if(val = rb_hash_aref(opts, ID2SYM(rb_intern("mu"))) != Qnil)mu = NUM2INT(val);

    iwork[0] = mu;
    iwork[1] = ml;
    lrw = max(20+16*neq, 22+10*neq+(2*ml+mu)*neq);
  }

  rwork = (double *)ALLOC_N(double, lrw);
  rwork[4] = h0;
  rwork[5] = hmax;
  rwork[6] = hmin;

  iwork[4] = ixpr;
  iwork[5] = max_step;
  iwork[6] = max_hnil;
  iwork[7] = max_ordn;
  iwork[8] = max_ords;

  // run lsoda
  lsoda_(ode_function, &neq, y, &t, &tout, &itol, &rtol, &atol, &itask, &istate, &iopt, rwork, &lrw, iwork, &liw, ode_jacobian_function, &jt);

  // check error
  if(istate<0){
    rb_raise(rb_eTypeError, "Something went wrong.");
  }

  for(i=0; i<neq; i++){
    VALUE num = DBL2NUM(y[i]);
    rb_ary_store(ret_arr, i, num);
  }

  xfree(y);
  xfree(iwork);
  xfree(rwork);
  return ret_arr;
}