archive/books/schelter/octave.bk

   1 \x06\x01\x19\x16\x05
   2 ((face (octave-eval 950 1106 1109 1343 1350 1506 1512 1675 1678 1956)) (book-command-arg))\f
   3             Sample Calls to Octave
   4
   5    If y is a function of t and satisfies the following differential
   6 equation
   7      y'' + .1* y' + sin(y) = 0;
   8
   9 If we use x(1) to denote a variable 'x' with subscript 1,
  10 then this may be coded in octave as a first order matrix equation:
  11
  12      y  = x(1),
  13      y' = x(2)
  14
  15 so eliminating y altogether:
  16
  17     x(1)' = x(2)
  18     x(2)' = -sin(x(1) - .1 * x(2)
  19
  20 In Octave we define a vector function 'pend' which encodes this, and
  21 then call the 'lsode' linear solver.  The t is ranging from values 0
  22 to 40, with 200 steps.  Then we plot the solution vector [x(1),x(2)]
  23 as two plots against t.  The [0.1,0.2] correspond to an initial value
  24 of X.  The initial value is at the first value of t (in this case 0)
  25 of the range requested.
  26
  27 The plot will have 2 curves in this case, and the 'line 1' corresponds
  28 to x(1), the 'line 2' corresponds to x(2).   Remember x = [y,y'] so
  29 that 'line 1' correpsonds to y and line 2 to y'.
  30
  31   function xdot = pend(x,t)
  32   xdot(1) = x(2); xdot(2) = -sin( x(1)) - 0.1*x(2); end
  33   sol=lsode( "pend",[0.1, 0.2], t = linspace(0,40, 200));
  34   plot( t, sol);
  35
  36
  37 function xdot = bruss(x,t)
  38 xdot(1) = 0;
  39 xdot(2) = 0;
  40 xdot(3) = x(1) - (x(2) + 1) * x(3) + x(3)*x(3)*x(4);
  41 xdot(4) = x(2) - x(3)*x(3)*x(4);
  42 endfunction
  43   sol=lsode( "bruss",[0.5, 1.5,0.5,1.5], t = linspace(0,40, 200));
  44   plot( t, sol);
  45
  46
  47
  48
  49   function xdot = pend(x,t)
  50   xdot(1) = x(2); xdot(2) = -sin( x(1)) - 0.1*x(2); end
  51   sol=lsode( "pend",[0.1, 0.2], t = linspace(0,20, 200));
  52   plot( t, sol);
  53
  54
  55
  56   function xdot = pend(x,t)
  57   xdot(1) = x(2);
  58   xdot(2) = -sin( x(1)) - 0.4*x(2)**2;
  59   end
  60   sol=lsode( "pend",[0.1, 0.2], t = linspace(0,20, 200));
  61   plot( t, sol);
  62
  63
  64 function xdot = pend(x,t)
  65   xdot(1) = .1*x(1) +.2* x(2)+.2*x(3);
  66   xdot(2) = .1*x(1) + 2*x(2) + .3*x(3);
  67   xdot(3) = .22*x(1)+x(2) + x(3);
  68   end
  69   sol=lsode( "pend",[1.0, 2.0,3.0], t = linspace(0.0,3.0, 200));
  70 /*  plot( t, sol(:,1))  plot just solution 1 */
  71   plot( t, sol);
  72
  73
  74 ;
  75