Search code examples
octavedifferential-equations

Differential equation; function in Octave


I try to solve a differential equation in Octave. In a first attempt, all my independent variables are set constant (n,y,w). This is my code:

function xdot= f(x,t);

  % xdot=zeros(1,1);

  X=1.44221E+12;
  IO=5.318E+11;
  GO=6.81E+11;
  beta=0;
  gamma=0;
  y=58.5021088;
  w= 31.29;
  n=1363.5;
  tw=0.4;
  tp=0.3;
  sw=0.07;
  sp=0.25;
  mw=0.593941689;
  mp=0.593941689;

  % aw=(1-tw)(sw+mw)
  % ap=(1-tp)(sp+mp) 

  xdot=-(X+IO*(1+gamma*(diff(n)/n+diff(y)/y))+GO*beta(tw(diff(n)/n+diff(w)/w)+tp(diff(y)/y-diff(w)/w)-(x*n)/((y-w)((1-tp)(sp+mp)+tp)+((1-tw)(sw+mw)+tw)*w))*x)/(IO*gamma+GO*beta*tw);


 endfunction

When I add

 t = [0:(1/360):10]

 x = lsode ("f", 39290000, t);

to solve the equation in the command line, I get this error:

error: index (0.843942): subscripts must be either integers 1 to (2^31)-1 or logicals
error: lsode: evaluation of user-supplied function failed
error: called from
f at line 23 column 7

It seems to me, in some way I misunderstood how to make the function. Any help?

EDIT:

This is my new code:

function xdot= f(x,t);

 X=1.44221E+12;
 IO=5.318E+11;
 GO=6.81E+11;
 beta=0;
 gamma=0;
 y=58.5021088;
 w= 31.29;
 n=1363.5;
 tw=0.4;
 tp=0.3;
 sw=0.07;
 sp=0.25;
 mw=0.593941689;
 mp=0.593941689;

 xdot=-X+IO-(x*n)/((y-w)*((1-tp)*(sp+mp)+tp)+w*(tw+(1-tp)*(sp+mp)))

endfunction

It does work if I copy it into the command line but if I start it as programm (f.m) and solve it then I get this error:

error: 'x' undefined near line 25 column 15
error: called from
f at line 25 column 7

Solution

  • From the documentation of lsode:

    Syntax: [x, istate, msg] = lsode (fcn, x_0, t)
    ...

    The first argument, fcn, is a string, inline, or function handle that names the function f to call to compute the vector of right hand sides for the set of equations. The function must have the form

    xdot = f (x, t)
    

    in which xdot and x are vectors and t is a scalar.

    So firstly, your f function is of the wrong form to the one expected by lsode, since it takes 5 arguments as opposed to 2.

    Secondly, inside your f function, the variables y, w, and n are overwritten, regardless of whether they were supplied, whereas t and x need to be supplied by the user.

    When lsode calls the handle f, it calls it with two arguments as per its specification, which means it instantiates t and y, and calls the function with the remaining input arguments undefined. Of the remaining arguments, w and n are given values inside the function, as we pointed out above, but x remains undefined. So when you run the final calculation which relies on x, you get an error saying it has not yet been defined during the call to this function.

    So, to fix this issue, rewrite your function such that its signature is f(x, t), which you can do so easily since y, w, and n are defined inside your function anyway, so they don't really need to be given as inputs.