Representation and Simulation of Dynamic Systems
- To: mathgroup at smc.vnet.net
- Subject: [mg56947] Representation and Simulation of Dynamic Systems
- From: "Caffa Vittorio Dr." <Caffa at iabg.de>
- Date: Wed, 11 May 2005 05:24:02 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
The behavior of (time-continuous, non-linear) dynamic systems can be
numerically investigated with NDSolve. One can first sketch a block
diagram of the system and then convert it into equations. Here is a toy
example after the conversion:
pos'[t] = vel[t]
vel'[t] = -k pos[t] + force[t] / m
This works fine if the variables are all states, as in the example
above. But often, in order to describe a given dynamic system you want
or you have to introduce some auxiliary variables (i.e. variables which
are not states). This is in fact the case if you want to describe a
generic dynamic system. Here are the standard equations:
x'[t] = f[x[t], u[t], t] (state equations)
y[t] = g[x[t], u[t], t] (output equations)
where: x = state vector, u = input vector, y = output vector, t = time.
In this case the components of the output vector are the "auxiliary"
variables.
I'm considering here a scheme for representing dynamic systems (possibly
using a block diagram as a starting point) which allows the usage of
auxiliary variables. This representation can be transformed into
equations for NDSolve automatically. After having solved the equations
it is possible to inspect not only the state variables but also the
auxiliary variables.
Comments or alternative solutions to the problem I'm considering are
welcome!
Procedure
o) Sketch the system on a piece of paper. Here is a toy example:
---------- [ -k ] ---------
| |
V |
force[t] --> [ 1/m ] --> + --> [ 1/s ] ---> [ 1/s ] ---> pos[t]
| |
| --------------> vel[t]
|
---------------------------> acc[t]
Note: [ 1/s ] is an integrator block
[ k ] is a gain block
o) Convert the sketch into a system description:
In[1]:= sys = {pos'[t] -> vel[t],
vel'[t] -> acc[t],
acc[t] -> -k pos[t] + force[t] / m};
Note: the arrow points to the source of the signal.
o) Make a list of the state variables:
In[2]:= states = {pos[t], vel[t]};
o) Form the differential equations (the following steps could be
performed by a function):
In[3]:= lhs = D[states, t]
Out[3]= {pos'[t], vel'[t]}
In[4]:= rhs = D[states, t] //. sys
force[t]
Out[4]= {vel[t], -------- - k pos[t]}
m
In[5]:= eqns = Join[Thread[lhs == rhs], {pos[0] == pos0, vel[0] ==
vel0}]
force[t]
Out[5]= {pos'[t] == vel[t], vel'[t] == -------- - k pos[t], pos[0] ==
pos0, vel[0] == vel0}
m
o) Specify the parameters:
In[6]:= params = {m -> 10, k -> 2, pos0 -> 0, vel0 -> 0, force[t] ->
Sin[t]};
o) Solve the differential equations:
In[7]:= sol = First[NDSolve[eqns /. params, states, {t, 0, 10}]]
Out[7]= {pos[t] -> InterpolatingFunction[{{0., 10.}}, <>][t],
vel[t] -> InterpolatingFunction[{{0., 10.}}, <>][t]}
o) Inspect the results (including auxiliary variables)
In[8]:= Plot[pos[t] /. sol, {t, 0, 10}]
Out[8]= -Graphics-
In[9]:= Plot[acc[t] //. sys /. params /. sol, {t, 0, 10}]
Out[9]= -Graphics-
Cheers, Vittorio
--------------------------------------------
Dr.-Ing. Vittorio G. Caffa
IABG mbH
Abt. VG 32
Einsteinstr. 20
85521 Ottobrunn / Germany
Tel. (089) 6088 2054
Fax: (089) 6088 3990
E-mail: caffa at iabg.de
Website : www.iabg.de
--------------------------------------------