/tags/cvs_final/octave-forge/extra/civil/inst/nlnewmark.m
MATLAB | 188 lines | 179 code | 8 blank | 1 comment | 6 complexity | 8c55ef118d9876d794f1f2dff90133fd MD5 | raw file
Possible License(s): GPL-2.0, BSD-3-Clause, LGPL-2.1, GPL-3.0, LGPL-3.0
- ## Copyright (C) 2000 Matthew W. Roberts. All rights reserved.
- ##
- ## This file is part of Octave.
- ##
- ## Octave is free software; you can redistribute it and/or modify it
- ## under the terms of the GNU General Public License as published by the
- ## Free Software Foundation; either version 2, or (at your option) any
- ## later version.
- ##
- ## Octave is distributed in the hope that it will be useful, but WITHOUT
- ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
- ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
- ## for more details.
- ##
- ## You should have received a copy of the GNU General Public License
- ## along with Octave; see the file COPYING. If not, write to the Free
- ## Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA 02111 USA.
- ## -*- texinfo -*-
- ## @deftypefn {Function File} {@var{x} =} nlnewmark(@var{m}, @var{c}, @var{q}, @var{f}, @var{dt}, @var{x0} = 0, @var{x'0} = 0, @var{alpha} = 1/2, @var{beta} = 1/4, @var{flags} = "")
- ## Computes the solution of non-linear second-order differential equations of the form
- ##
- ## @example
- ## @var{m} @var{x}'' + @var{c} @var{x}' + @var{q}(@var{x}, @var{x}', @var{x}'') = @var{f}
- ## @end example
- ##
- ## where @var{x}' denotes the first time derivative of @var{x} and @var{q} is
- ## a non-linear function.
- ##
- ## If the function is called without the assigning a return value
- ## then @var{x} is plotted versus time.
- ##
- ## @strong{Inputs}
- ##
- ## @table @var
- ## @item m
- ## The mass of the body.
- ## @item c
- ## Viscous damping of the system.
- ## @item q
- ## The name of a function
- ## that returns the value of the resisting force for a given
- ## displacement. The form of Q must be:
- ##
- ## @example
- ## @var{F} = Q( @var{x})
- ## @end example
- ##
- ## where @var{F} is the restoring force for the state vector @var{x}
- ## = [@var{u}, @var{u'}, @var{u''}]; displacement
- ## (@var{u}), velocity (@var{u'}), and acceleration.
- ## @item f
- ## The forcing function as a time sampled or impulse vector
- ## (see @strong{Special Cases}).
- ## @item dt
- ## The time step -- assumed to be constant
- ## @item x0
- ## Initial displacement, default is zero
- ## @item x'0
- ## Initial velocity, default is zero
- ## @item alpha
- ## Alpha Coefficient -- Controls "artificial damping" of the system.
- ## Unless you have a really good reason, this should be 1/2 which is
- ## the default.
- ## @item beta
- ## Beta Coefficient -- This coefficient is used to estimate the form of the
- ## system acceleration between time steps. Values between 1/4 and 1/6 are
- ## common. The default is 1/4 which is unconditionally stable.
- ## @item flags
- ## A string value which defines special cases. The cases are defined by
- ## unique characters as explained in @strong{Special Cases} below.
- ## @end table
- ##
- ## @strong{Outputs}
- ##
- ## @table @var
- ## @item x
- ## Matrix of size (3, @code{length(@var{f})}) with time series of displacement
- ## (@var{x}(1,:)), velocity (@var{x}(2,:)), and acceleration (@var{x}(3,:))
- ## @end table
- ##
- ## @strong{Special Cases}
- ##
- ## The @var{flags} variable is used to define special cases of analysis as
- ## follows.
- ##
- ## "i" - Impulse forcing function. The forcing function, @var{f} is a
- ## vector of impulses instead of a sampled time history.
- ## "n" - The stiffness is non-linear. In this case, @var{k} is a string
- ## which contains the name of a function defining the non-linear
- ## stiffness.
- ##
- ## @end deftypefn
- ## Author: Matthew W. Roberts
- ## Created: May, 2000
- function x = nlnewmark( M, C, Q, f, dt, u0, v0, alpha, beta, flags)
- global nlnewmark_status;
- # Take care of unititalized input variables
- if( nargin < 6)
- u0 = 0;
- endif
- if( nargin < 7)
- v0 = 0;
- endif
- if( nargin < 8)
- alpha = 0.5;
- endif
- if( nargin < 9)
- beta = 0.25;
- endif
- if( nargin < 10)
- flags = "";
- endif
- # check for flags
- if( findstr( flags, "i"))
- _local_impulse = 1; # local variable
- else
- _local_impulse = 0;
- endif
- # xxBEGINxx
- # initialize x
- x = zeros( 3, length(f));
- x(1,1) = u0;
- x(2,1) = v0;
- # compute the initial acceleration
- if( _local_impulse)
- # an initial impulse has the effect of instantaneously changing
- # the initial velocity.
- v0 = v0 + f(1)/M;
- x(2,1) = v0;
- # Now, initial acceleration comes from solving m*a0 + c*v0 + Q = 0
- ## :FIXME: here I'm assuming that the acceleration is not used in Q
- x(3,1) = ( - C * v0 - feval( Q, [u0, v0, 0]) ) / M;
- else
- # Compute a0 from m*a0 + c*v0 + Q = f0
- ## :FIXME: here I'm assuming that the acceleration is not used in Q
- x(3,1) = ( f(1) - C * v0 - feval( Q, [u0, v0, 0]) ) / M;
- endif
- # define some constants so that we won't need to recalculate each loop:
- c1 = dt^2 * (0.5 - beta);
- c2 = dt * (1 - alpha);
- # These need to be known by __nlnewmark_fcn__
- nlnewmark_status.a1 = dt^2*beta;
- nlnewmark_status.a2 = dt*alpha;
- nlnewmark_status.C = C;
- nlnewmark_status.M = M;
- nlnewmark_status.Q = Q;
- rhs(3) = 0; # default value - always this for impulse
- for i = 2:length(f)
- % create the rhs
- nlnewmark_status.rhs(1) = x(1, i-1) + dt * x(2, i-1) + c1 * x(3, i-1);
- nlnewmark_status.rhs(2) = x(2, i-1) + c2 * x(3, i-1);
- if( ! _local_impulse)
- nlnewmark_status.rhs(3) = f(i);
- endif
- # solve for x - using previous value as starting guess...
- [x(:, i), info] = fsolve( "__nlnewmark_fcn__",
- [x(1, i-1); x(2, i-1); x(3, i-1)]);
- info
- # add the impulse effect...
- if( _local_impulse)
- x(2,i) = x(2,i) + f(i)/M;
- endif
-
- endfor
- if( nargout < 1)
- t = 0:dt:(length(f)-1)*dt;
- length(t);
- plot( t, x(1,:));
- endif
- endfunction