B. Background of the invention Models We consider simulation of systems that consist of a group of subsystems. Especially we consider multiphysics systems; term multiphysics system refers to a system whose subsystems may belong to different physical domains. A subsystem can be, for example, a body, an actuator, an electronic component, or a hydraulic valve. The subsystems are connected to each other with algebraic constraints that can, for example, enforce the continuity of a hydraulic flow, or specify the relative position or orientation of two rigid bodies. Under suitable assumptions these constraints can be taken into account using Lagrange multipliers [Asch98, Schwe99]. Often, a constraint involves only few of the state variables of a submodel. For a description on multiphysics modeling, see [Lay98].

Mathematically such systems can be described as ordinary differential equations with constraints, which is called a differential-algebraic equation (DAE). The most common form of DAEs is F (x' (t), x (t), t) =O, where the matrix Fx (,), subscript denoting a partial derivative and prime denoting a time derivative, is singular. DAEs arising from the physical domains usually have a more distinct structure. For example, DAE models of articulated rigid multibody dynamics with constraints are of the semi-explicit form y'(t) = I (y (t), z (t)) z'(t)=k(y(t),z(t),v(t))

0=p (y (t)), (1) where the vectors y (t) and z (t) refer to the positions and velocities of the bodies, and v (t) is a vector of Lagrange multipliers related to the position constraints p (y (t)) = 0, see [Fet80, Hai96, Asch98, Sha98].

Multiphysics models may also contain velocity constraints, that is, the algebraic constraints are of the form p (y (t), z (t)) = 0.

We assume these constraints are componentwise holonomic or Pfaffian. By holonomic we refer to constraints involving position variables only, and constraints involving velocity variables that can be integrated to obtain a constraint involving positions only. Pfaffian constraints refer to such conditions imposed on infinitesimal displacements that are transformable into integrable velocity variable constraints.

The assumption is needed to justify the use of Lagrange multipliers in taking the constraints into account. For details see, e. g., [Lay98]. Inequality constraints, also referred to as unilateral constraints, see [Pfe96, Bro99], can be handled separately, if needed.

A differential algebraic equation can be classified according to its index. For different index definitions, see [Hai96]. We mention here the differentiation index that is defined in the literature as the minimal number of differentiations with respect to time that is needed for both sides of a DAE to obtain a differential equation for all the variables. For example, the DAE (1) is of index 3, since three differentiations in time are needed to obtain a differential equation for v (t) as well.

A standard technique in the numerical solution of a DAE is to reduce its index by a suitable technique. One of the most important index reduction schemes is the Gear- Gupta-Leimkuhler (GGL) formulation [Gea85], that can be used to transform an index-3 problem into an index-2 problem. Semi-explicit index-3 DAE reduced into index-2 form using GGL reduction is of the form x'(t)=k(x(t),u(t),t) g (x (t), t) =D,

where x(t)=[yT(t),zT(t)]T and u (t) = [vT (t), nT (t)] T and n (t) is the vector of Lagrange multipliers resulting from the GGL index reduction. For the precise functional form offand g see, e. g. , [Hai96].

In the following, we focus on semi-explicit index-2 problems. Problems of lower differentiation index can be regarded as special cases of these, and index-3 problems can be reduced to index-2 problems. Problems of practical significance with an order higher than 3 seldom appear in multiphysics. Most efficient numerical DAE solvers, in terms of complexity and computing time, are designed for index-2 problems. For certain higher-index DAEs some methods, though quite cumbersome, exist; see, e. g. , [Wen98].

Differential-algebraic equation solvers Under certain conditions the algebraic constraints could be eliminated by moving to generalized coordinate representation (e. g. [Schw99] ), but the process is complicated and results in a very complex set of differential equations for the remaining coordinates. In practice, solution methods are based on solving the model above.

Two different approaches, direct discretization and stabilization, can be distinguished. Direct discretization approaches [Hai96], either one-step or multistep, are implicit integration schemes. They discretize the dynamics of the problem in time and solve a system of nonlinear equations, specified by (2), for the values of x (t) and u (t) at each timestep. Newton's method is usually applied. Semi-implicit schemes replace the solution of the nonlinear equation by a solution of a set of suitably selected linear equations.

In stabilization methods, the DAE is first transformed into an ordinary differential equation with an invariant manifold by reducing its differentiability index to zero.

This is done by differentiating the algebraic constraints until a differential equation is obtained for the algebraic variables, too. The ordinary differential equation is then solved with any suitable numerical integration scheme, and stabilized onto the invariant manifold. Baumgarte's stabilization [Bau72] attempts to stabilize the original continuous-time equation. Recent approaches [Asch97, Asch98] stabilize the discretized system at the end of each integration timestep by a projection onto the manifold. Even for direct discretization methods it may be preferable to lower the index of the DAE by suitable differentiations and modifications.

Real-time simulation and its requirements Here, simulation of a physical system refers to solving the DAE model above repeatedly in time. Simulation can be performed as a batch process, if all the inputs of the model are known beforehand. In real-time simulation, it is required that the DAE model advances synchronously with the global time. Since a real-time simulation responds, by definition, to external events or data within a predefined deadline, it allows, for example, control system testing and tuning and operational training in a practically riskless manner. Real-time simulation of physical systems is an inherent part of modern design, as it allows virtual prototyping at all the stages of the design cycle without the need for expensive actual prototypes.

Real-time simulation is feasible, if the computing time required to numerically solve the differential equation models for a time step is shorter than the time step itself. Typically, the simulator needs time for the data visualization and communication as well. Hence, in real-time simulation, numerical solution of the differential equations can take only a fraction of the time step. As the number of subsystems in the physical system increases, the dimension of the differential- algebraic equation becomes larger. Since the computation time is always proportional to the dimension of the system (typically, to its 3rd power), the real time property of the simulation is eventually lost at some point, no matter how fast a processor is being used.

Parallel and distributed simulation In attempts to speed up computation, one can either increase the speed of a single processor, or increase the number of processors in the underlying computer system and to parallelize the computations by distributing the computational workload among them in a suitable way. The former approach is always limited by the current technological limits. In the latter approach, the decrease in the computational time, or the speedup, can in principle be made arbitrarily large. This is one reason why parallel and distributed computing today forms an active field of research. Other sources of motivation for distributed computing, given in [FujOO], are the possibility to geographical distribution, possibility to connect different computer architectures without program code porting, and increased fault tolerance.

An important concept in measuring the performance of a system operating in parallel is the speedup S. It is defined as the time it takes a program to execute in serial divided by the time it takes to execute in parallel. In practice, the achievable speedup depends on the design of the numerical algorithm and its applicability to parallelization. Amdahl's law (e. g. [Ber89] ) states that the largest possible speedup S depends on the algorithm as follows: S= (s +p)/(s +p/N).

Here, N is the number of processors, s is the time a serial processor spends on the part of the code that cannot be parallelized, and p is the time a serial processor spends on the code part that can be parallelized. It is implicitly assumed thatp does not depend on N. The law shows that it is the algorithm and not the number of processors that limits the speedup.

We call an algorithm scalable if the speedup predicted by Amdahl's law is significant. Otherwise we call the algorithm inscalable. In computer science, scalability is often defined as the ability of a computer application or product to continue to function well as the application or product itself or its context is changed in size or volume. In practice, this means that if more than one processor becomes available, the computational workload of the application can be distributed among the units. Clearly, a scalable application requires a scalable numerical algorithm.

Significance of communication overhead Usually the part of the code that is parallelized must communicate with other parts of the code. When parallel program segments are executed in separate processors with separate memories, which is the case in distributed computing, like in PC clusters, the time needed in the communication becomes dominant. It is therefore important that different segments of parallel applications communicate as little as possible.

Inscalability of current DAE solvers Present numerical solution algorithms of DAEs are centralized, designed for one processor, and hence inherently inscalable. The algorithms do allow some speedup, e. g. , through parallelization of elementary matrix operations, but the algorithms per

se remain essentially sequential. There does not seem to exist a scalable DAE solver.

C. Summary of the invention An object of the invention is to present a method for parallel and distributed simulation of a physical system consisting of subsystems that are connected to each other with constraints that are taken into account using Lagrange multipliers.

The object is achieved with a two-level algorithm that decouples the integration of the subsystems, that is, the determination y (t) and z (t), and the determination of the Lagrange multipliers u (t).

A method according to the invention is characterized by that which is specified in the characterizing part of an independent claim directed to a method for simulating a physical system consisting of a plurality of subsystems. A simulation arrangement according to the invention is characterized by that which is specified in the characterizing part of an independent claim directed to a simulation arrangement, and computing units according to the invention are characterized by that which is specified in the characterizing parts of independent claims directed to computing units. A computer program set according to the invention and computer programs according to the invention are characterized by that which is specified in the characterizing parts of independent claims directed to a computer program set and to computer programs.

The dependent claims describe some preferred embodiments of the invention.

At a time step, each subsystem is first integrated separately at the lower level from the given current state. According to the invention, this integration can be performed for each subsystem independently of other subsystems. At the upper level, the Lagrange multipliers are then adjusted such that the constraints connecting the subsystems stay satisfied at the end of the time step. If needed, the adjusted values of the multipliers are updated to the lower level, and the integration of the subsystems is repeated. In a simulation method in accordance with the invention, subsystem integration and Lagrange multiplier value selection is typically carried out in each simulation time step. Without loss of generality we speak of integration even though a subsystem may have only algebraic dependencies.

Parallelism is achieved through the fact that with given, constant values of the Lagrange multipliers, and linear connecting constraints, the subsystems can be integrated independently of each other. For nonlinear constraints the same holds when the connecting states of the neighboring systems are extrapolated suitably. Hence, the integration of the subsystems can be distributed to any number of processors less than or equal to the amount of subsystems.

Each subsystem can be integrated with an integration method that is most suitable for it. As mentioned above, the subsystem can also be an algebraic equation or even a DAE. Hence, recursive organization of the computation is also possible. The only requirement is that the subsystems must be connected to each other with constraints which can be taken into account using Lagrange multipliers. Holonomic or Pfaffian constraints are examples of such constraints.

At the upper level, the Lagrange multipliers that are communicated to the lower level are selected by solving an equation (typically nonlinear) with a suitable numerical method. The dimension of the equation is typically much smaller than that of the DAE. Thus, significant scalability should be achieved.

Speedup of the computation is, however, only one of the benefits of the method. In a decentralized environment, also the construction and preparation of the simulation model can be distributed, which leads to conceptually new ways to organize simulations. A large number of individuals or other actors, like computers, can be granted access to a distributed simulation. With a commonly agreed interface, the participants can insert or remove models in the simulation and observe the behavior of the overall system. A significant application area for such distributed simulations is virtual environments (see [FujOO]).

In the simulation, the subsystems preferably only exchange connecting states and Lagrange multipliers. In other words, a subsystem does not have to reveal its internal structure. Given the initial state of the system, the Lagrange multiplier values corresponding to the constraints between subsystems, and a timestep length at some time instant, a subsystem only has to announce its connecting state at the end of the timestep. This means, for example, that competing product suppliers can share a common simulation without exposing the details of the product implementation. In a larger perspective, a common simulation interface allows rapid distributed simulation-based prototyping. For example, potential buyers could

compare the products of different suppliers by constructing simulation models over a network. Since the suppliers can also provide the computational capacity for the simulation of their subsystem, the buyer does not necessarily need particularly excessive computational power.

To summarize, according to the invention - the integration or, more generally, time-stepping of the subsystems is decomposed from the determination of the Lagrange multipliers and from each other - the Lagrange multipliers are determined by solving a system of equations by any suitable method.

D. Brief description of figures The invention is below described in more detail with reference to the following figures.

Figure 1 presents the decomposition principle and communication between the subsystems and the multiplier coordinator, Figure 2 presents the physical system treated in the numerical example, Figure 3 presents a stroboscopic image of the motion of the example system, Figure 4 presents the computed Lagrange multiplier trajectories and compares them with a reference solution, Figure 5 presents the state variable trajectories computed by the method together with a reference solution, Figure 6 presents the error in the algebraic constraints in the numerical example, and Figure 7 presents schematically a simulation arrangement according to the invention.

E. Determination of the multipliers

As defined earlier, simulation in the context of the invention refers to solving a differential-algebraic equation, that in general is of the form (2). Assume the system to be simulated consists of m subsystems, each described as X, (t) f hi(xi(t))=0, i=1,...,m, where Sri is a Lagrange multiplier vector related to the internal constraints hi (dim hi=nhi) of subsystem i. In general, the constraints h ; may be absent for some or all subsystems. Next, let us divide the states of subsystem i into two groups, the internal states xiI(t) and the connecting states xiC(t), and organize the state vector such that xi(t)=[xiI(t),xiC(t)]. Typically, the number of the connecting states, nic, is much smaller than the number of the internal states. Consider then vector-valued, independent constraints that connect each subsystem to its neighbors as gk (xC (t)) =O, k=l,..., p. where xC (t) = [x1C(t),...,xmC(t)]. Denote the dimension of gk (. ) by ngk. To avoid overconstrainedness, we require that the total number of constraints, ng = ng,+... +ngp, be strictly less than the number of degrees of freedom of the unconnected subsystems. Then, the model of the whole system can be expressed as Here xi(t) # Rni, #i(t) # Rnhi, u(t) = [u1(t),...,up(t)] # Rng, and g (x (t)) = [g1T(xc(t),...,gpT(xc(t))]T. The simulation is started at to from a given consistent initial state x- [x, T(t0),...,xmT(t0)]T, v0=[v1T(t0),...,vmT(t0)]T, Uo= [u1T(t0),...,upT(t0)]T, and continues over a time interval [totl], t, < An initial state is consistent if it satisfies the algebraic constraints and the constraints of the form

where I is one less than the differentiation index of the DAE. For details on consistent initial states, see e. g. [Schw99].

Let Fi be a DAE solver or, in the absence of the constraints/ ! an ODE (ordinary differential equation) integrator for the subsystem i. For notational convenience we assume that Fi is a one-step method and the integration timestep At is constant. Also multistep methods and variable stepsize within a solver can be used. Then, at time j#t, xi ( (j + 1) At) = Fi (xi (j#t),#Ci(xj#t)),#Ci(u(j#t))), where I-Ic (x (t)) contains the connecting states of the subsystems that are connected to subsystem i, and #Ci (u (t)) contains the multipliers corresponding to the connecting constraints of the subsystem i. As Figure 1 illustrates, u (j#t) and x c (iAt) are communicated in step 102 for time-stepping the subsystems i. During At in step 103 in Figure 1, u (j#t) is held constant. If the computation has been successful, we have g (x (j#t) = [g1(xC(j#t))T,...,gk(xC(j#t))T]T = 0. However, this does not automatically (even theoretically) hold for g (xC((j+1)#t)). In the equation above, xi(j#t) and rIc (x (t)) are considered as a fixed initial condition and, besides At and Fi, x ((j+l) t) depends only on u(j#t). Thus, we can select u (j#t) such that g (x ((j+1)#t))=0. This selection is done in step 101 in Figure 1, and from the subsystem time-stepping the subsystem-specific connecting states are communicated for this selection in step 104. In other words, we solve a nonlinear equation <BR> <BR> <BR> <BR> <BR> Xj ((j+1)#t)=Fi(xi(j#t),#Ci(x(j#t)),#Ci(u(j#t))), i=1,...,m, (4)<BR> <BR> <BR> g(xC((j+1)#t)=0.

The solution methods for nonlinear equations are iterative, which means that the methods approach the solution stepwise using knowledge on the previous step to take the next step. Perhaps the most well known iterative method for nonlinear equations is the Newton iteration and its variants. Also other solution methods can be used. When Newton iteration is applied, the time stepping procedure for time step j can be described as follows.

1. Fix an initial estimate u° (jet). Use e. g. the value from the previous timestep or a suitably extrapolated value. Set k:=0.

2. Time-step subsystems with Fi using uk (j#t) to obtain xk ((j + 1) #t).

3. If llg (XC, k ((j+1)#t))# is suitably small, accept uk(j#t) and xk ((j+1)#t).

Proceed to the next timepoint.

Otherwise, take a Newton step; uk+1(j#t)=uk(j#t)-[gu(xC((j+1)#t))]-1 g(xC((j+1)#t)).

Set k:=k+1 and return to 2.

Newton iteration The matrix go (.) can be computed analytically using the chain rule of differentiation if the partial derivatives Fi,u are known. Nevertheless, in distributed computation, this is not necessarily a realistic assumption. In that case, a numerical scheme must be applied. In Newton iteration, The Jacobian is estimated anew at each step k as Here zlui is a vector of zeros with ith element replaced by a small perturbation and u is an abbreviation for u(j#t). The fixed states xC (j#t) have been omitted, and g (xC ((j+1)#t)) is considered only as a function of the multipliers u. Note that to obtain the approximations, xc ((j+1)#t) has to be evaluated for each perturbed u.

In simplified Newton iteration, the Jacobian of step k is reused at next steps, and it is updated only every nth step. The simplification reduces the average computational load but is likely to increase the number of steps due to worse convergence.

We finally mention the possibility to avoid even the computation of the finite differences by the use of secant updates to approximate the Jacobian. Letting the superscript k denote the Newton iteration step number at timestep i, an iteration with the Broyden rank-1 update, given e. g. in [Sto83], reads dk:=(Bk)-1g(uk) <BR> <BR> uk+1 :-uk _kdk<BR> <BR> <BR> <BR> <BR> <BR> <BR> pk :=uk+1-uk,qk:=g(uk+1)-g(uk)<BR> B pk q P p<BR> <BR> <BR> <BR> <BR> <BR> Bk+1:=Bk###(qk-Bkpk)pkT.<BR> <BR> <P> T<BR> <BR> <BR> p@p It can be seen that after a suitable initialization, the Jacobian approximation B can be computed using previous iterates and function values at these iterates only. The stepsize A can be computed, e. g. , by minimizing the norm of g at each step. It should be noted that the special structure of the problem probably permits the derivation of other, more specific iteration methods as well.

F. Integration of the subsystems As seen above, the differential equations in (3) can be integrated separately with a given u (t) using any suitable integration method. To show how the integration can be distributed, we use the assumption that the system consists of loosely connected subsystems that exchange information only via the connecting constraints and corresponding Lagrange multipliers.

Consider first the case when the last term in the 1st row of equation (3) does not contain xC, jxi. A sufficient condition for this is that the constraints gk are linear in x. Then, given the multiplier vector u, and the initial state, subsystem i can be time- stepped independently of the other subsystems.

When the constraints gk are nonlinear, connecting states of subsystem j can be present in the model of subsystem i. The separation then requires that the

connecting states of subsystem j are extrapolated. For short enough time intervals, linear extrapolation, i. e., xc (t+dt) =xc. (t) +xc' (t) dt, dt<At can be used. Note that xj' (t) need not be computed separately, as it is readily available in the right hand side of the differential equation in (3). Also quadratic extrapolation can be used, but it requires an additional computational effort.

G. Distributed model management Model management here refers to the actions that affect the structure of the simulation model. Examples of such actions are submodel insertion and removal, or a change in the submodel's connecting structure. It also includes organizing a global access to the submodels to produce, e. g. , a snapshot on the system's state at any time instant, and the actions with which the upper-level coordinator conducts the simulation. Such actions include issues like stepsize selection and integration accuracy requirement.

Unsupervised and distributed model management requires that there exists a commonly agreed interface, or a description on how submodels communicate with the coordinator and vice versa. Despite the middleware used to implement the system, at least the following items typically have to be included in the interfaces: 1. Coordinator interface A reference to the coordinator itself Notification mechanism through which a submodel announces its appearance or disappearance in the model Initial connecting state information of the neighboring states at each time step Extrapolation information of the neighboring connecting states at each time step Current stepsize Accuracy requirements.

2. Subsystem time stepper interface Connecting constraint computation Connecting state output

Whole subsystem state output for a snapshot Local time stepping error estimate Participation in the computation of consistent initial values In the case that models include unilateral constraints, a mechanism to inform a constraint violation.

H. Communication capacity requirements Let us assume that a given system consists of m subsystems and, on the average, 1. each subsystem has n state variables, of which nc are connecting variables 2. Each subsystem is connected to r neighboring subsystems with k linear constraints in each connection 3. The method used in solving the Lagrange multipliers converges in q iterations 4. Computing one timestep takes At, which is much larger than the lags related to the communication channel, and much less than the simulation timestep At.

Assume, furthermore, that a globally consistent initial state is provided to each subsystem.

In the simulation at the time jZJt, j E N, the upper level coordinator must provide each subsystem with the rk multipliers, which means transmitting rkm floating point numbers (fpn). The subsystems return a tentative connecting state (nCm fpn), after which the coordinator communicates the updated multiplier values to the subsystems. The cycle occurs q times, until the subsystems return the state that <BR> <BR> <BR> satisfies g (. ) =0 with a pre-specified accuracy. Thus, with these assumptions, the total number of floating point numbers transmitted per one second is A = q (rkm + nCm) qm (rk + nC) At At Taking q=2, r=2, nC= 4 and k=3, taking m=10 and At=0. 03 seconds gives the capacity requirement of approximately 6700 fpn/s. Assuming 32 bit arithmetic, this corresponds to slightly more than 200 kbit/s for the overall traffic. For an individual integration unit, the capacity requirement is 20 kbit/s, since m does not contribute in it. Such a capacity can be achieved using an ordinary telephone line modem.

It can be seen that the more connected the system is in terms of the average number of connecting constraints k, their average dimension r, and the average dimension of the subsystems'connecting variables n, the more communication is needed.

Increasing the timestep decreases the capacity requirement. It should be noted, however, that q is likely to increase with increasing At. Note that A does not depend on n.

If all the constraints are nonlinear, each subsystem needs to know the connecting states of the neighboring r systems, and their time derivatives for the 1st order extrapolation purposes. This means the transmission of additional 2rnCm fpn at each timestep. Taking the same numerical values as above, we end up with the overall capacity requirement of 550 kbit/s, or 55 kbit/s for a single time stepping or integration unit.

I. Speedup estimates The simplest integration scheme for DAEs, implicit Euler method (e. g. [Asch98]) solves a system of nonlinear equations at each timestep. The unknowns in the equation are the state variables and the Lagrange multipliers. Using the nomenclature above, the dimension of the equation system is thus approximately (nm+mrkl2), the two appearing from the fact that each constraint involves two bodies. One step of Newton iteration involves solving a system of linear equations.

When dense matrix algebra is used, it takes a time proportional to the 3rd power of the dimension of the system. Thus, the time required for a single Newton step is hence proportional to (nm+mrkl2) 3, provided that an analytical expression for the Jacobians of the state equation and the constraints are available. Numerical approximation of the Jacobians takes a time proportional to the dimension of the problem, and the time required for a single Newton step becomes proportional to nm+mrkl2 + (nm+mrkl2) 3.

In the method described here, the subsystems are integrated either simultaneously or in parallel, the multipliers are iterated, and the process is repeated until convergence is achieved. The nonlinear equation to be solved involves only the Lagrange multipliers and has the dimension mrkl2. At each Newton step, the subsystems are integrated. When explicit integration is used, this takes a time proportional to nm, the total number of state variables. If, furthermore, the Jacobian is approximated with finite differences at each step, each of the m subsystems has to be integrated

additional rk times, which takes a time proportional no nmrk. Hence, one step with the method takes a time proportional to nm+nmrk+ (mrkl2) 3.

The parallelizable part of the code is the integration of the subsystems. Maximal speedup is obviously achieved when the number of processors is m+1 (one for the Newton iteration and m for the integration). Thus, the speedup for one Newton iteration is <BR> <BR> <BR> S=nm+mrk/2+(nm+mrk/2)³ (1+2n)3.<BR> <BR> <BR> <BR> n + nrk + (mrk/2)³ m-># rk Taking n=13, m=10 and the same parameter values as above gives an estimate S = 152. If we assume that an equal number of Newton steps is required in both methods, this is the speedup of the integration as well.

In this comparison, the speedup is always greater than one. Increasing connectivity increases the number of Lagrange multipliers and reduces the speedup, whereas increasing average subsystem state dimension increases it.

J. Example We consider a double pendulum with a spring attached to the ceiling and the lower mass, see figure 2. The system is in a gravity field and the masses are subject to simple viscous friction. Subsystems are modeled using Lagrangian DAEs as follows: Pendulum i, i=1, 2 In this Section J, x and y are the coordinates of the masses mi, c denotes the viscous friction coefficient and g is the acceleration due to gravity.

Spring

For convenience we assume there is a (small) mass m at the end of the spring; x and y denote its position and Fx, y the components of the spring force F, where k is the spring constant and rl is the rest length.

The subsystems are connected to each other with the constraints x12+y12-l12=0 (x2-x- (Y2 2 X3-X2 = 0 y3-yw=0.

Here//=1, 2 are the lengths of the pendulums. The interpretation of the constraints is that the first pendulum is attached to the origin, the second pendulum is attached to the first one, and the spring is attached to the second pendulum. We use the parameter values given in Table 1 below.

Simulation results We simulate the system from the initial state x (0) =0.76, yl (0) =1. 19, x2 (0) =x3 (0) =1. 52, y2 (0) =y3 (0) =0, from 0 to 5 seconds. The subsystems are integrated with ODE23 of Matlab [Anon01] using relative precision of 10-7 and absolute precision of 10-8. A reference solution is computed by eliminating the Lagrange multipliers corresponding to the constraints and using Baumgarte stabilization [Bau72] with a real double root at - 10. For illustration, a stroboscope picture of the motion of the masses is presented in Figure 3. MI 20 kg m2 10 k m3 1 kg 11 1. 4142 m 12 1. 4142 m c 3 N/m/s g 9.81 . mls2 k 30 N/m rl 2 m (xf,yf) (5,0) Table 1. The parameter values

The trajectories of the Lagrange multipliers uj (t), j=1... 4, corresponding to the constraints above, are presented in Figure 4 with At=0. 02 seconds. The reference solution is presented in Figure 4 with smooth lines. The subsystems are integrated separately, and a linear extrapolation of the states of the adjacent subsystems is used in the nonlinear blocks of G. The extrapolation does not overly affect the accuracy.

The state variable trajectories are shown in Figure 5. They coincide with the reference solution, which is also presented in Figure 5, within the drawing accuracy.

The 2-norm of the constraint error g is plotted in Figure 6. When u ((t) is used as the initial estimate for u (jAt) in the iteration, only two Newton steps for each i are needed to achieve the accuracy of 10-5. Tests with simplified Newton iteration show that the number of iterations increase to 3-4. Simplified Newton iteration would be a preferable choice since computing the Jacobian by finite differences takes most of the computation time.

At this stage, the computations are carried out with a single processor, and hence judgments on the computing times are not speculable. We estimate that when implemented properly, this distributed simulation method should be applicable also to real-time simulation.

K. Simulation arrangement and simulation software

Figure 7 illustrates schematically a simulation arrangement 700 according to the invention. The simulation arrangement 700 comprises a first computing unit 701, which is responsible for Lagrange multiplier value selection, and at least one second computing unit 702 for time-stepping a subsystem to a tentative next subsystem state. A subsystem state is formed of values of subsystem-specific variables and a subsystem connecting state is formed of values of subsystem-specific connecting variables, which are involved in the constraints.

The first computing unit 701 comprises - means for communicating with said at least second computing unit 702, and - means for selecting values for Lagrange multipliers, said means arranged to accept tentative next subsystem states corresponding to the tentative next subsystem connecting states as next subsystem states, if the constraints are valid for the tentative subsystem connecting states and arranged otherwise to select new values for Lagrange multipliers for sending to the at least one second computing unit.

A second computing unit 702 comprises - means for communicating with said first computing unit 701, - means for storing a current subsystem state, and - means for time-stepping a subsystem using Lagrange multiplier values received from the first computing unit and starting from a current subsystem state and current subsystem connecting state (s) of neighboring subsystem (s) connected to a specific subsystem.

The first and second computing units may be, for example, PC computers or other computers connected to each other using any communication network and communicating with each other using any standard or proprietary communication protocol. Alternatively, they may be processors in a multiprocessor computer and communicate with each other using means and protocols specific to the multiprocessor computer.

Preferably the first and second computing units are realized by using suitable software in combination with conventional computing devices.

It is possible that a second computing unit is responsible for time-stepping more than one subsystem of the simulated physical system. Furthermore, it is possible that a first computing unit and a second computing unit are realized in one device (computer or processor). It is furthermore possible that the second computing units

exchange information about the connecting states of the subsystems between themselves, but as the first computing unit also needs the information about the connecting states, it is typically advisable to transmit this information between the second computing units via the first computing unit.

In view of the foregoing description, it will be evident to a person skilled in the art that various modifications may be made within the scope of the invention. While some preferred embodiments of the invention have been described in detail, it should be apparent that many modifications and variations thereto are possible, all of which fall within the scope of the invention as defined by the appended independent claims.

References [AnonO1] Matlab Reference Manual, Mathworks Inc. , 2001.

[Asch97] Ascher, U. ,"Stabilization of invariants of discretized differential systems", Numerical Algorithms Vol. 14 pp. l-24, 1997.

[Asch98] Ascher, U. and Petzold, L., Computer Methods for Differential Equations and Differential-Algebraic Equations, Society for Industrial and Applied Mathematics, Philadelphia, 1998.

[Bau72] Baumgarte, J. ,"Stabilization of constraints and integrals of motion in dynamical systems", Computer Methods in Applied Mechanics, Vol. 1, pp. 1-16, 1972.

[Ber89] Bertsekas, D. and Tsitsiklis, J. , Parallel and Distributed Computing : Numerical Methods, Prentice-Hall, Englewood Cliffs, 1989.

[Bro99] Brogliato, B., Nonsmooth Mechanics, 2"d edition, Springer, London, 1999.

[Fet80] Fetter, A. and Walecka, J., Theoretical Mechanics of Particles and Continua, McGraw-Hill, 1980.

[FujOO] Fujimoto, R., Parallel and Distributed Simulation Systems, Wiley, New York, 2000.

[Gea85] Gear, C. , Gupta, G. , and Leimkuhler, B. ,"Automatic Integration of Euler- Lagrange Equations with Constraints", Journal of Computational and Applied Mathematics, Vol. 12-13,1985, pp. 77-90.

[Hai96] Hairer, E. and Wenner, G., Solving Ordinary Differential Equations II- Stiff and Differential-Algebraic Problems, 2nd ed. , Springer, Berlin, 1996.

[Lay98] Layton, R., Principles of Analytical System Dynamics, Springer, New York, 1998.

[Pfe96] Pfeiffer, F. and Glocker, C., Multibody dynamics with unilateral contacts, Wiley, New york, 1996.

[Schw99] Schwerin, R., Multibody System Simulation, Springer, Berlin, 1999.

[Sha98] Shabana, A., Dynamics of Multibody Systems, Cambridge Unversity Press, Cambridge, 1998.

[Sto80] Stoer, J. and Bulirsch, R., Introduction to Numerical Analysis, Springer, New York, 1980.

[Wen98] Wensch, J. ,"An eight stage fourth order partitioned Rosenbrock method for multibody systems in index-3 formulation", Applied Numerical Mathematics Vol. 27 No. 2 pp. 171-183,1998.