octave_packages/odepkg-0.8.2/odepkg_examples_ode.m

   1 %# Copyright (C) 2008-2012, Thomas Treichl <treichl@users.sourceforge.net>
   2 %# OdePkg - A package for solving ordinary differential equations and more
   3 %#
   4 %# This program is free software; you can redistribute it and/or modify
   5 %# it under the terms of the GNU General Public License as published by
   6 %# the Free Software Foundation; either version 2 of the License, or
   7 %# (at your option) any later version.
   8 %#
   9 %# This program is distributed in the hope that it will be useful,
  10 %# but WITHOUT ANY WARRANTY; without even the implied warranty of
  11 %# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12 %# GNU General Public License for more details.
  13 %#
  14 %# You should have received a copy of the GNU General Public License
  15 %# along with this program; If not, see <http://www.gnu.org/licenses/>.
  16
  17 %# -*- texinfo -*-
  18 %# @deftypefn {Function File} {[@var{}] =} odepkg_examples_ode (@var{})
  19 %# Open the ODE examples menu and allow the user to select a demo that will be evaluated.
  20 %# @end deftypefn
  21
  22 function [] = odepkg_examples_ode ()
  23
  24   vode = 1; while (vode > 0)
  25     clc;
  26     fprintf (1, ...
  27       ['ODE examples menu:\n', ...
  28        '==================\n', ...
  29        '\n', ...
  30        '   (1) Solve a non-stiff "Van der Pol" example with solver "ode78"\n', ...
  31        '   (2) Solve a "Van der Pol" example backward with solver "ode23"\n', ...
  32        '   (3) Solve a "Pendulous" example with solver "ode45"\n', ...
  33        '   (4) Solve the "Lorenz attractor" with solver "ode54"\n', ...
  34        '   (5) Solve the "Roessler equation" with solver "ode78"\n', ...
  35        '\n', ...
  36        '   Note: There are further ODE examples available with the OdePkg\n', ...
  37        '         testsuite functions.\n', ...
  38        '\n', ...
  39        '   If you have another interesting ODE example that you would like\n', ...
  40        '   to share then please modify this file, create a patch and send\n', ...
  41        '   your patch with your added example to the OdePkg developer team.\n', ...
  42        '\n' ]);
  43     vode = input ('Please choose a number from above or press <Enter> to return: ');
  44     clc; if (vode > 0 && vode < 6)
  45       %# We can't use the function 'demo' directly here because it does
  46       %# not allow to run other functions within a demo.
  47       vexa = example (mfilename (), vode);
  48       disp (vexa); eval (vexa);
  49       input ('Press <Enter> to continue: ');
  50     end %# if (vode > 0)
  51   end %# while (vode > 0)
  52
  53 %!demo
  54 %! # In this example the non-stiff "Van der Pol" equation (mu = 1) is
  55 %! # solved and the results are displayed in a figure while solving.
  56 %! # Read about the Van der Pol oscillator at
  57 %! #   http://en.wikipedia.org/wiki/Van_der_Pol_oscillator.
  58 %!
  59 %! function [vyd] = fvanderpol (vt, vy, varargin)
  60 %!   mu = varargin{1};
  61 %!   vyd = [vy(2); mu * (1 - vy(1)^2) * vy(2) - vy(1)];
  62 %! endfunction
  63 %!
  64 %! vopt = odeset ('RelTol', 1e-8);
  65 %! ode78 (@fvanderpol, [0 20], [2 0], vopt, 1);
  66
  67 %!demo
  68 %! # In this example the non-stiff "Van der Pol" equation (mu = 1) is
  69 %! # solved in forward and backward direction and the results are
  70 %! # displayed in a figure after solving. Read about the Van der Pol
  71 %! # oscillator at http://en.wikipedia.org/wiki/Van_der_Pol_oscillator.
  72 %!
  73 %! function [ydot] = fpol (vt, vy, varargin)
  74 %!   ydot = [vy(2); (1 - vy(1)^2) * vy(2) - vy(1)];
  75 %! endfunction
  76 %!
  77 %! vopt = odeset ('NormControl', 'on');
  78 %! vsol = ode23 (@fpol, [0, 20], [2, 0], vopt);
  79 %! subplot (2, 3, 1); plot (vsol.x, vsol.y);
  80 %! vsol = ode23 (@fpol, [0:0.1:20], [2, 0], vopt);
  81 %! subplot (2, 3, 2); plot (vsol.x, vsol.y);
  82 %! vsol = ode23 (@fpol, [-20, 20], [-1.1222e-3, -0.2305e-3], vopt);
  83 %! subplot (2, 3, 3); plot (vsol.x, vsol.y);
  84 %!
  85 %! vopt = odeset ('NormControl', 'on');
  86 %! vsol = ode23 (@fpol, [0:-0.1:-20], [2, 0], vopt);
  87 %! subplot (2, 3, 4); plot (vsol.x, vsol.y);
  88 %! vsol = ode23 (@fpol, [0, -20], [2, 0], vopt);
  89 %! subplot (2, 3, 5); plot (vsol.x, vsol.y);
  90 %! vsol = ode23 (@fpol, [20:-0.1:-20], [-2.0080, 0.0462], vopt);
  91 %! subplot (2, 3, 6); plot (vsol.x, vsol.y);
  92
  93 %!demo
  94 %! # In this example a simple "pendulum with damping" is solved and the
  95 %! # results are displayed in a figure while solving. Read about the
  96 %! # pendulum with damping at
  97 %! #   http://en.wikipedia.org/wiki/Pendulum
  98 %!
  99 %! function [vyd] = fpendulum (vt, vy)
 100 %!   m = 1;    %# The pendulum mass in kg
 101 %!   g = 9.81; %# The gravity in m/s^2
 102 %!   l = 1;    %# The pendulum length in m
 103 %!   b = 0.7;  %# The damping factor in kgm^2/s
 104 %!   vyd = [vy(2,1); ...
 105 %!          1 / (1/3 * m * l^2) * (-b * vy(2,1) - m * g * l/2 * sin (vy(1,1)))];
 106 %! endfunction
 107 %!
 108 %! vopt = odeset ('RelTol', 1e-3, 'OutputFcn', @odeplot);
 109 %! ode45 (@fpendulum, [0 5], [30*pi/180, 0], vopt);
 110
 111 %!demo
 112 %! # In this example the "Lorenz attractor" implementation is solved
 113 %! # and the results are plot in a figure after solving. Read about
 114 %! # the Lorenz attractor at
 115 %! #   http://en.wikipedia.org/wiki/Lorenz_equation
 116 %! #
 117 %! # The upper left subfigure shows the three results of the integration
 118 %! # over time. The upper right subfigure shows the force f in a two
 119 %! # dimensional (x,y) plane as well as the lower left subfigure shows
 120 %! # the force in the (y,z) plane. The three dimensional force is plot
 121 %! # in the lower right subfigure.
 122 %!
 123 %! function [vyd] = florenz (vt, vy)
 124 %!   vyd = [10 * (vy(2) - vy(1));
 125 %!          vy(1) * (28 - vy(3));
 126 %!          vy(1) * vy(2) - 8/3 * vy(3)];
 127 %! endfunction
 128 %!
 129 %! A = odeset ('InitialStep', 1e-3, 'MaxStep', 1e-1);
 130 %! [t, y] = ode54 (@florenz, [0 25], [3 15 1], A);
 131 %!
 132 %! subplot (2, 2, 1); grid ('on');
 133 %!   plot (t, y(:,1), '-b', t, y(:,2), '-g', t, y(:,3), '-r');
 134 %!   legend ('f_x(t)', 'f_y(t)', 'f_z(t)');
 135 %! subplot (2, 2, 2); grid ('on');
 136 %!   plot (y(:,1), y(:,2), '-b');
 137 %!   legend ('f_{xyz}(x, y)');
 138 %! subplot (2, 2, 3); grid ('on');
 139 %!   plot (y(:,2), y(:,3), '-b');
 140 %!   legend ('f_{xyz}(y, z)');
 141 %! subplot (2, 2, 4); grid ('on');
 142 %!   plot3 (y(:,1), y(:,2), y(:,3), '-b');
 143 %!   legend ('f_{xyz}(x, y, z)');
 144
 145 %!demo
 146 %! # In this example the "Roessler attractor" implementation is solved
 147 %! # and the results are plot in a figure after solving. Read about
 148 %! # the Roessler attractor at
 149 %! #   http://en.wikipedia.org/wiki/R%C3%B6ssler_attractor
 150 %! #
 151 %! # The upper left subfigure shows the three results of the integration
 152 %! # over time. The upper right subfigure shows the force f in a two
 153 %! # dimensional (x,y) plane as well as the lower left subfigure shows
 154 %! # the force in the (y,z) plane. The three dimensional force is plot
 155 %! # in the lower right subfigure.
 156 %!
 157 %! function [vyd] = froessler (vt, vx)
 158 %!   vyd = [- ( vx(2) + vx(3) );
 159 %!          vx(1) + 0.2 * vx(2);
 160 %!          0.2 + vx(1) * vx(3) - 5.7 * vx(3)];
 161 %! endfunction
 162 %!
 163 %! A = odeset ('MaxStep', 1e-1);
 164 %! [t, y] = ode78 (@froessler, [0 70], [0.1 0.3 0.1], A);
 165 %!
 166 %! subplot (2, 2, 1); grid ('on');
 167 %!   plot (t, y(:,1), '-b;f_x(t);', t, y(:,2), '-g;f_y(t);', \
 168 %!         t, y(:,3), '-r;f_z(t);');
 169 %! subplot (2, 2, 2); grid ('on');
 170 %!   plot (y(:,1), y(:,2), '-b;f_{xyz}(x, y);');
 171 %! subplot (2, 2, 3); grid ('on');
 172 %!   plot (y(:,2), y(:,3), '-b;f_{xyz}(y, z);');
 173 %! subplot (2, 2, 4); grid ('on');
 174 %!   plot3 (y(:,1), y(:,2), y(:,3), '-b;f_{xyz}(x, y, z);');
 175
 176 %# Local Variables: ***
 177 %# mode: octave ***
 178 %# End: ***