[Octave-cvsupdate] SF.net SVN: octave:[6196] trunk/octave-forge/main/control/inst/initial.m From: - 2009-09-01 18:58 ```Revision: 6196 http://octave.svn.sourceforge.net/octave/?rev=6196&view=rev Author: paramaniac Date: 2009-09-01 18:58:04 +0000 (Tue, 01 Sep 2009) Log Message: ----------- initial: support for discrete systems Modified Paths: -------------- trunk/octave-forge/main/control/inst/initial.m Modified: trunk/octave-forge/main/control/inst/initial.m =================================================================== --- trunk/octave-forge/main/control/inst/initial.m 2009-09-01 17:17:50 UTC (rev 6195) +++ trunk/octave-forge/main/control/inst/initial.m 2009-09-01 18:58:04 UTC (rev 6196) @@ -53,10 +53,10 @@ ## @seealso{impulse, lsim, step} ## @example ## @group -## . -## x = A x x(0) = x0 +## . +## Continuous Time: x = A x , y = C x , x(0) = x0 ## -## y = C x +## Discrete Time: x[k+1] = A x[k] , y[k] = C x[k] , x[0] = x0 ## @end group ## @end example ## @end deftypefn @@ -64,8 +64,7 @@ ## Author: Lukas Reichlin ## Created: August 16, 2009 ## based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel -## Version: 0.1 alpha -## TODO: Needs thorough testing for discrete systems! +## Version: 0.1 function [y_r, t_r, x_r] = initial (sys, x_0, t_final, dt) @@ -183,16 +182,29 @@ endfor if (nargout == 0) # plot information - for k = 1 : n_out - subplot (n_out, 1, k) - plot (t, y(:, k)) - grid on - if (k == 1) - title ("Response to Initial Conditions") - endif - ylabel (sprintf ("Amplitude %s", sysgetsignals (sys, "out", k, 1))) - endfor - xlabel ("Time [s]") + if (digital) # discrete system + for k = 1 : n_out + subplot (n_out, 1, k) + stairs (t, y(:, k)) + grid on + if (k == 1) + title ("Response to Initial Conditions") + endif + ylabel (sprintf ("Amplitude %s", sysgetsignals (sys, "out", k, 1))) + endfor + xlabel ("Time [s]") + else # continuous system + for k = 1 : n_out + subplot (n_out, 1, k) + plot (t, y(:, k)) + grid on + if (k == 1) + title ("Response to Initial Conditions") + endif + ylabel (sprintf ("Amplitude %s", sysgetsignals (sys, "out", k, 1))) + endfor + xlabel ("Time [s]") + endif else # return values y_r = y; t_r = t; This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6203] trunk/octave-forge/main/control/inst/initial.m From: - 2009-09-02 14:10 ```Revision: 6203 http://octave.svn.sourceforge.net/octave/?rev=6203&view=rev Author: paramaniac Date: 2009-09-02 14:10:11 +0000 (Wed, 02 Sep 2009) Log Message: ----------- initial: fix bug, add test Modified Paths: -------------- trunk/octave-forge/main/control/inst/initial.m Modified: trunk/octave-forge/main/control/inst/initial.m =================================================================== --- trunk/octave-forge/main/control/inst/initial.m 2009-09-02 13:35:05 UTC (rev 6202) +++ trunk/octave-forge/main/control/inst/initial.m 2009-09-02 14:10:11 UTC (rev 6203) @@ -64,7 +64,7 @@ ## Author: Lukas Reichlin ## Created: August 16, 2009 ## based on __stepimp__.m of Kai P. Mueller and A. Scottedward Hodel -## Version: 0.1 +## Version: 0.2 function [y_r, t_r, x_r] = initial (sys, x_0, t_final, dt) @@ -96,7 +96,7 @@ endif ## code adapted from __stepimp__.m - if (nargin < 3) + if (nargin < 4) ## we have to compute the time when the system reaches steady state ## and the step size eigw = eig (sys2ss (sys)); @@ -125,28 +125,34 @@ endfor if (nk == 0) - t_final = 10; + if (nargin < 3) + t_final = 10; + endif dt = t_final / 1000; else eigw = eigw(find (eigw)); eigw_max = max (abs (eigw)); dt = 0.2 * pi / eigw_max; - eigw_min = min (abs (real (eigw))); - t_final = 5.0 / eigw_min; + if (nargin < 3) + eigw_min = min (abs (real (eigw))); + t_final = 5.0 / eigw_min; - ## round up - yy = 10^(ceil (log10 (t_final)) - 1); - t_final = yy * ceil (t_final / yy); + ## round up + yy = 10^(ceil (log10 (t_final)) - 1); + t_final = yy * ceil (t_final / yy); + endif - n = t_final / dt; + if (! digital) + n = t_final / dt; - if (n < 50) - dt = t_final / 50; - endif + if (n < 50) + dt = t_final / 50; + endif - if (n > 2000) - dt = t_final / 2000; + if (n > 2000) + dt = t_final / 2000; + endif endif endif endif @@ -169,7 +175,7 @@ ## preallocate memory y = zeros (l_t, n_out); - x_vec = zeros (l_t, n_st); + x_arr = zeros (l_t, n_st); ## make sure that x is a row vector x = reshape (x_0, length (x_0), 1); @@ -177,7 +183,7 @@ ## simulation for k = 1 : l_t y(k, :) = C * x; - x_vec(k, :) = x; + x_arr(k, :) = x; x = F * x; endfor @@ -208,11 +214,68 @@ else # return values y_r = y; t_r = t; - x_r = x_vec; + x_r = x_arr; endif endfunction -## TODO: Add a purely continuous system test -## TODO: Add a purely discrete system test +%!shared initial_c, initial_c_exp, initial_d, initial_d_exp +%! +%! A = [ -2.8 2.0 -1.8 +%! -2.4 -2.0 0.8 +%! 1.1 1.7 -1.0 ]; +%! +%! B = [ -0.8 0.5 0 +%! 0 0.7 2.3 +%! -0.3 -0.1 0.5 ]; +%! +%! C = [ -0.1 0 -0.3 +%! 0.9 0.5 1.2 +%! 0.1 -0.1 1.9 ]; +%! +%! D = [ -0.5 0 0 +%! 0.1 0 0.3 +%! -0.8 0 0 ]; +%! +%! x_0 = [1, 2, 3]; +%! +%! sysc = ss (A, B, C, D); +%! sysd = c2d (sysc, 2); +%! +%! [yc, tc, xc] = initial (sysc, x_0, 0.2, 0.1); +%! initial_c = round (100 * [yc, tc, xc]) / 100; +%! +%! [yd, td, xd] = initial (sysd, x_0, 4); +%! initial_d = round (100 * [yd, td, xd]) / 100; +%! +%! yc_exp = [ -1.00 5.50 5.60 +%! -0.99 5.09 5.77 +%! -0.95 4.69 5.76 ]; +%! +%! tc_exp = [ 0.00 +%! 0.10 +%! 0.20 ]; +%! +%! xc_exp = [ 1.00 2.00 3.00 +%! 0.59 1.69 3.09 +%! 0.24 1.52 3.10 ]; +%! +%! initial_c_exp = [yc_exp, tc_exp, xc_exp]; +%! +%! yd_exp = [ -1.00 5.50 5.60 +%! -0.65 3.17 4.22 +%! -0.54 2.62 3.50 ]; +%! +%! td_exp = [ 0 +%! 2 +%! 4 ]; +%! +%! xd_exp = [ 1.00 2.00 3.00 +%! -0.42 1.52 2.32 +%! -0.35 1.25 1.93 ]; +%! +%! initial_d_exp = [yd_exp, td_exp, xd_exp]; +%! +%!assert (initial_c, initial_c_exp) +%!assert (initial_d, initial_d_exp) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```
 [Octave-cvsupdate] SF.net SVN: octave:[6221] trunk/octave-forge/main/control/inst/initial.m From: - 2009-09-06 08:08 ```Revision: 6221 http://octave.svn.sourceforge.net/octave/?rev=6221&view=rev Author: paramaniac Date: 2009-09-06 08:08:29 +0000 (Sun, 06 Sep 2009) Log Message: ----------- initial: rework test Modified Paths: -------------- trunk/octave-forge/main/control/inst/initial.m Modified: trunk/octave-forge/main/control/inst/initial.m =================================================================== --- trunk/octave-forge/main/control/inst/initial.m 2009-09-05 22:46:14 UTC (rev 6220) +++ trunk/octave-forge/main/control/inst/initial.m 2009-09-06 08:08:29 UTC (rev 6221) @@ -244,36 +244,38 @@ %! sysd = c2d (sysc, 2); %! %! [yc, tc, xc] = initial (sysc, x_0, 0.2, 0.1); -%! initial_c = round (100 * [yc, tc, xc]) / 100; +%! initial_c = round (1e4 * [yc, tc, xc]) / 1e4; %! %! [yd, td, xd] = initial (sysd, x_0, 4); -%! initial_d = round (100 * [yd, td, xd]) / 100; +%! initial_d = round (1e4 * [yd, td, xd]) / 1e4; %! -%! yc_exp = [ -1.00 5.50 5.60 -%! -0.99 5.09 5.77 -%! -0.95 4.69 5.76 ]; +%! ## expected values computed by the "dark side" %! -%! tc_exp = [ 0.00 -%! 0.10 -%! 0.20 ]; +%! yc_exp = [ -1.0000 5.5000 5.6000 +%! -0.9872 5.0898 5.7671 +%! -0.9536 4.6931 5.7598 ]; %! -%! xc_exp = [ 1.00 2.00 3.00 -%! 0.59 1.69 3.09 -%! 0.24 1.52 3.10 ]; +%! tc_exp = [ 0.0000 +%! 0.1000 +%! 0.2000 ]; %! +%! xc_exp = [ 1.0000 2.0000 3.0000 +%! 0.5937 1.6879 3.0929 +%! 0.2390 1.5187 3.0988 ]; +%! %! initial_c_exp = [yc_exp, tc_exp, xc_exp]; %! -%! yd_exp = [ -1.00 5.50 5.60 -%! -0.65 3.17 4.22 -%! -0.54 2.62 3.50 ]; +%! yd_exp = [ -1.0000 5.5000 5.6000 +%! -0.6550 3.1673 4.2228 +%! -0.5421 2.6186 3.4968 ]; %! %! td_exp = [ 0 %! 2 %! 4 ]; %! -%! xd_exp = [ 1.00 2.00 3.00 -%! -0.42 1.52 2.32 -%! -0.35 1.25 1.93 ]; +%! xd_exp = [ 1.0000 2.0000 3.0000 +%! -0.4247 1.5194 2.3249 +%! -0.3538 1.2540 1.9250 ]; %! %! initial_d_exp = [yd_exp, td_exp, xd_exp]; %! This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```