 # Dynamic Deformation of Uniform Elastic Two-Layer Objects

This thesis presents a two-layer uniform facet elastic object for real-time simulation based on physics modeling method. It describes the elastic object procedural modeling algorithm with particle system from the simplest one-dimensional object, to more complex two-dimensional and three-dimensional objects. The double-layered elastic object consists of inner and outer elastic mass spring surfaces and compressible internal pressure. The density of the inner layer can be set different from the density of the outer layer; the motion of the inner layer can be opposite to the motion of the outer layer. These special features, which cannot be achieved by a single layered object, result in improved imitation of a soft body, such as tissue's liquidity non-uniform deformation. The construction of the double-layered elastic object is closer to the real tissue's physical structure. The inertial behavior of the elastic object is well illustrated in environments with gravity and collisions with walls, ceiling, and floor. The collision detection is defined by elastic collision penalty method and the motion of the object is guided by the Ordinary Differential Equation computation. Users can interact with the modeled objects, deform them, and observe the response to their action in real time.

## Authors

##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

### 1.1 Definitions

#### 1.1.1 Deformable Object

In engineering mechanics, “deformable object” refers to an object whose shape can be changed due to an applied force, such as tensile (pulling), compressive (pushing), bending, or tearing forces. The deformation can be categorized as the following, depending on the types of material and the forces applied:

• Elastic deformation (small deformation) is reversible. The object shape is temporarily deformed when tension is applied and it returns to its original shape when force is removed. An object made of rubber has a large elastic deformation range; silk cloth material has a moderate elastic deformation range; crystal has almost no elastic deformation range.

• Plastic deformation (moderate deformation) is not reversible. The object shape is deformed when tension is applied and its shape is partially returned to its original form when the force is removed. Objects such as silver and gum, which can be stretched at their original length, cannot completely restore their original shapes after deformation.

• Fracture deformation (large deformation) is not reversible, but is different from the plastic deformation. The object is permanently deformed when it is irreversibly bent, torn, or broken apart after the material has reached the end of the elastic deformation ranges. All materials will experience fracture deformation when sufficient force is applied.

#### 1.1.2 Elastic Object

Elastic objects belong to a subset of soft body deformable objects. They are dynamic objects that change shape significantly and keep constant volume in response to collision. They can be bent, stretched, and squeezed. Moreover, they restore their previous shape after deformation. Elastic objects can be divided into two domains:

• Large elastic deformation, such as fluid deformation, which focuses on flows through space. It tracks velocity and material properties at fixed points in space.

• Small elastic deformation, such as tissue deformation, which uses particle systems to identify chunks of matter and track their position, acceleration, and velocity.

Within this wide research range of soft body simulation, this work has focused on small elastic deformable object simulation, such as tissue animation. Even though there has been many valuable contribution related to this field, there are still many difficulties in accomplishing to realistic and efficient deformable simulation.

### 1.2 Animation Techniques

This section introduces some basic concepts related to the elastic simulation, such as the subject animation method. Animation relies on persistence of vision and refers to a series motion illusions resulting from the display of static images in rapid-shown succession. In traditional animation, squash and stretch are exaggerated for elastic objects. In order to be efficient when working with many of single frame images (or simply frames), inbetweening and cel animation [TJ84] have been introduced by Disney for manual traditional animation.

The rate of the animation refers to how many frames are displayed within a given amount of time. If the rate is too low, which is lower than the brain visual retention, the animation becomes jerky because the brain retains the empty frame from the previous image to the next image.

A frame rate, which is the time between two updates of the display, describes the update frequency. In computer games, frames are often discussed in terms of frames per second (fps). The lower bound for smooth animation is between 22 to 30 frames per second.

For many years’ research, computer-animation has been developed dramatically to replace the amount of manual traditional animation. The techniques of key-framing, morphing, and motion capture [HO99] have been widely used.

• Key-frame animation: is based on manual animation. It is a sequence of images of the same object with its local transformations, e.g. values for translation, rotation and scale, computed by interpolating between keyframes.

• Morphing: is a method usually used to estimate and generate a sequence of frames between one object to the other object with same number of vertices. Morphing is a good animation technique when using skeletal animation would be too complex, e.g. facial animation.

• Motion capture: is the method that attaches sensors on actors bodies and records the data for their movements and apply these data to a computer generated characters.

### 1.3 Elastic Animation

There are two different methods about elastic animation modeling, which depends on the predefined simulation or simulation in real time.

###### Kinematic modeling

predefines the positions and velocities of objects. It does not concern what causes movement and how things get where they are in the first place and only deals with the actual movement. For example, given that a ball’s initial speed is 10 kilometers per hour on a perfect smooth plane, we can use kinematic method to calculate how far it travels in two hours.

###### Dynamic modeling

also termed as physically based modeling, is the study of masses and forces that cause the kinematic quantities, such as acceleration, velocity, and position, to change as time progresses. For example, when we know the ball’s initial speed, we need to know how far it travels after an external force dynamically applied to it.

For elastic object movement, the dynamic methods calculate how the soft body behaves after external force applied dynamically. The animator does not need to specify the exact path of an object compared to using the kinematic modeling method. The system predefines the initial condition of the elastic object, such as position and gravity force. The animation of the object movement is updated each time step based on the acceleration derived from Newton’s Laws of motion. The dynamic simulation method is more advanced, easier to achieve the realistic motion than kinematic method. Therefore, we will only represent dynamic simulation of elastic object in this thesis.

### 1.4 Applications

Elastic modeling has been developed and used in several fields, including geometric modeling, computer vision, computational mathematics, physics engines, bio-mechanics, engineering, character animation, and many other fields

[GM97].

###### Character Animation

There is much advanced animation modeling software, which has advanced features for modeling, texturing, and lighting. However, for modeling the simulation of elastic objects, 3D artists have to do it manually, frame-by-frame because most of the current 3D software does not provide soft object simulation functionality. Artists have to use not only their drawing skills and intuition, but also posses some knowledge of physics to make the objects behave as if they are in the real world or close.

The techniques of the non-physically based modeling of the elastic object include modeling the group of control points and changing their property parameters manually frame-by-frame. The virtual objects will not convince audiences because no natural laws of physics are applied. Moreover, key-frame animation is an inefficient way to model elastic objects without functionality provided by software. Hence, most of 3D film animators have to ignore the movement details of soft objects.

The latest version of the most advanced animation tool, Maya, provides the Soft Mod Deformer tool, which allows smooth sculpting of a group of objects [Wag07]. However, users need to have knowledge about how to use this complex software in order to access this advanced functionality. Moreover, users can only animate elastic object with Kinematic modeling method by setting values through the software interface rather than interact with the object in real time.

###### Computer Games

Compared to the fancy and lifelike character animation widely used in 3D films, computer games are more concerned about computation efficiency because users interact with the software in real time. As one might notice, the majority of computer games do not portray the characters in detail, even with the mesh and texture modeling. It is not likely that elastic simulation will be widely introduced to computer games because existing elastic models usually require expensive calculation and are inconvenient to use in real time simulation.

###### Surgical Training

Surgeons benefit from the rapid development of computer graphics and hardware techniques. The computer generated visual virtual environment imitates the reality of medical operations and organ construction to fulfill the training purpose. This new application improves surgical outcomes and decreases the research costs. However, the reality and accuracy of the software always require high-end knowledge of physics, mathematic and heavy computation. It makes it difficult for users to interact with virtual objects in real time.

### 1.5 Thesis Goal

The elastic object for dynamic simulation, which has been widely used, is the one layer elastic surface with different content within. The soft objects can be squashed and stretched according to external and internal forces applied on them. Its computation depends on geometric modeling methods and physical equations. However, this method is inefficient to imitate the behaviors of real human’s tissue because human’s or animal’s soft body does not consist of only one layer with either liquid or air inside. Figure 1.2 from a biological research group demonstrates the complexity of human tissue [McE05]. A tissue is composed of epidermis, dermis, fat, fascia, and muscle layer.

• The epidermis is skin’s outermost layer. It is responsible for the skin coloring because it contains the skin’s pigment.

• The dermis, which is the layer that lies below the epidermis, consists entirely of living cells. It provides the skin elasticity because this layer is composed of bundles of fibers and small blood vessels.

• The fat, the fascia, and the muscle layer are hypodermis layer of skin. It is an inner layer of and cushion for the skin. This subcutaneous tissue layer varies throughout the body region and hormonal influence. The fat and muscle increase the tension of the skin and protect the bones.

Soft tissues are multi-composite layers; therefore, one layer elastic object is not accurate to model the kind of soft body exemplified by human tissue. Moreover, it is difficult to represent the object’s inertia caused by the internal material realistically and its liquidity motion based on the various material densities.

In this thesis, we investigate a more accurate two-layer elastic object. The outer layer of the elastic object represents the epidermis and the dermis layer of a real tissue. The inner layer of this object corresponds to the hypodermis layers of skin. This two-layer computer generated elastic object provides more accurate modeling method based on the main feature of human tissue. Its deformation is based on the realistic physical consistency of tissues and the laws of established physics. The objective of this new model is to be convincing and to have distinct realism to the animated scene by applying proper physics. The program should be easy in implementation, convenient to re-use, and gives best elastic body behavior at the minimum cost rather than model the absolute complex object with the exact accurate physical equations. Users should be able to interact with the soft body in real-time and the collision detection and response must be handled correctly.

### 1.6 Organization

This chapter starts with the introduction of elastic objects, their applications, some basic concepts related to physical based deformable simulation, and the thesis goal. Chapter 2 surveys and analyzes the existing elastic simulation system and its problems. Chapter 3 describes the modeling methodology of elastic objects in one-dimension, two-dimension, and three-dimension. Physically-based modeling and simulation map a natural phenomena to a computer simulation program. There are two basic processes in this mapping: mathematical modeling and numerical solution [Lin06]. Chapter 4 introduces mathematical modeling, which describes natural phenomena by mathematical equations. Chapter 5 presents the dynamics numerical equation of motion by using ODE (ordinary differential equation) associated with our elastic simulation system, and discusses the complexity and improvement of the different numerical time integrator of Euler, Midpoint, and Runge Kutta 4th order. Chapter 6 presents the detailed design and implementation of the simulation system. Chapter 7 shows our experimental results with the animation sequences of the elastic object simulation and discusses the effects of changing the simulation parameters. Chapter 8 gives the conclusion of our current system, summarizes our contributions based on the existing elastic simulation models, and analyzes the possible development and related work in the future.

### 2.1 Existing Elastic Object Models

###### Particle Model

has been used by Reeves [Ree83] and to model the natural phenomena such as fire, water, liquid as shown in Figure 2.1. Particles move under the force field and constraint without interacting with each other.

The advantage of this particle model is that the method is easy to implement. It is the earliest model to imitate the natural phenomena.

The disadvantage is that all the particles are independent and there are no forces connecting the particles. Therefore, for some phenomena, such as springs and mass, the method cannot achieve.

###### Deformable Surface

was introduced first time by Terzopoulos et al. [TPBF87], using finite difference for the integration of energy-based Lagrange equations based on Hooke’s law.

It was very successful in creating and animating surfaces. However, this method requires not only the discretization of the surfaces into spline patches, but also the specification of local connectivity for spring mass systems. These involve a significant amount of manual preprocessing before the surface model can be used.

###### Linear Mass Spring System

has been widely used for modeling elastic objects as shown in Figure 2.2. It is actually derived from the particle model; however, it simplifies the modeling of the inter-particle connection by using flexible springs. Three dimensional systems contain a finite set of masses connected by springs, which are assumed to obey Hooke’s Law.

This method was first introduced by Terzopoulos [TJF89] to describe melting effects. Particles, which are connected by springs, have an associated temperature as one of their properties. The stiffness of the spring is dependent on the temperature of the linked two particles. Increased temperature decreases the spring stiffness. When the temperature reaches the melting point, the stiffness becomes zero.

The advantages of mass-spring model are that it is easy to construct and display the simulation dynamically.

The disadvantages are that such system restricts to only the objects with small elastic deformation with approximation of the true physics, not for the objects that require large elastic deformation, such as fluid. This method also has difficulties to simulate the separation and fusion of a constant volume object. Moreover, the spring stiffness is problematic. If the spring is too weak, for the closed shape model with only simple springs to model the surface will be very easy to collapse. If we try to avoid the collapse, we need to model with spring stiffer, and then we will have difficulty to choose the elasticity because the object is nearly rigid. Another disadvantage is that the mass spring system has less stability and requires the numerical integrator to take small time steps [DW92] than FEM model.

###### Finite Element Method

known as FEM Model [GM97]

, is the most accurate physical model compared to other models. It treats deformable object as a continuum, which means the solid bodies with mass and energy distributed all over the object. This continuum model is derived from equations of continuum mechanics. The whole model can be considered as the equilibrium of a general object subjected to external forces. The deformation of the elastic object is a function of these forces and the material property. The object will stop deformation and reach the equilibrium state when the potential energy is minimized. The applied forces must be converted to the associated force vectors and the mass and stiffness are computed by numerically integrating over the object at each time step, so the re-evaluation of the object deformation is necessary and requires heavy pre-processing time

[GM97].

The advantage of FEM model is that it gives more realistic deformation result than mass-spring system because the physics are more accurate.

The disadvantage is that the system lacks efficiency. Because the energy equation will be used, the FEM is only efficient for the small deformation of the elastic object, such as application to the plastic material, which has a small deformation range. Alternatively, the object has less control elements needed to be computed, as in cloth deformation. If we need to simulate the human soft body parts or facial animation, the deformation rate is very high. It will be very difficult and sometimes impossible to carry out the integration procedure over the entire body. Therefore, it has been limited to apply in real-time system because of the heavy computational effort (usually it is done off-line). Moreover, the implementation is complicated.

###### Fluid Based Model

[DL02] consists of two components: an elastic surface and a compressible fluid as shown in Figure 2.3. The surface is represented as a mass spring system. The fluid is modeled using finite difference approximations to the Navier-Stokes equations of fluid flow. Figure 2.3 describes how this model simulates the fluid flows down a sink simulated. The inner layer is modeled by a particle system, which is similar to real water molecules. Using the numerical methods, the motion of each particle can be computed. In this example, the motion of the each particle is at the center of the basin, and points down to the sink.

The fluid based model uses physically based modeling and it produces realistic fluid animation. It illustrates the behavior of fluid in environments with gravity and collisions with planes.

The disadvantage of this model is the heavy computation for both elastic surface and density inside fluid. It also provides a solution to the constant volume problem.

###### Pressure Model

was introduced by M. Matyka [Mat03, Mat04b, Mat04b]. It simulates an elastic deformable object with a internal pressure based on the ideal gas law as shown in Figure 2.4. The object volume is calculated approximately by bounding box, shaped as sphere, cube, or ellipsoid. Another method to determine the object volume is based on Gauss’s Theorem.

Advantage of this model is that it gives very convincing effects about elastic property in real time simulation. The object behaves like a balloon filled only with air.

However, it can not imitate more interesting effects, such as human tissue. It can not achieve the effect of semi-liquid deformable object because the air pressure density is uniform inside of the object, which is different from liquid with non-uniform density. It is not accurate for describing the inertia of the semi-liquid object.

### 2.2 Summary of The Existing Models

Previous work on deformable object animation uses physically-based methods with local and global deformations applied directly to the geometric models. Based on the survey of the existing elastic models, we conclude their usage as the two types:

• Interactive models are used when speed and low latency are most important and physical accuracy is secondary. Typical examples include mass-spring models and spline surfaces used as deformable models. These can satisfy the character animation with exaggerated unrealistic deformation.

• Accurate models are chosen when physical accuracy is important in order to accomplish the surgical training purpose which requires the accurate tissue modeling. The continuum simulation model, for instance, the most accurate model, FEM, is not ideal for simulation requiring real time interaction and the object undergoing large deformation.

In short, elastic object simulation is a dilemma of demanding accuracy and interactivity.

### 3.1 Graphics Objects Modeling Methods

###### Polygonal Methods

create geometric objects that can be described by their convex planar polygonal surfaces. These methods are easy to describe the shapes of mathematical objects rendered on graphics system. However, they are difficult to describe physical objects, such as cloud and fire, and their complex behaviors combined with physical laws [Ang03].

###### Procedural Methods

build natural phenomena, 3D models and textures in an algorithmic manner and generate polygons only during the rendering process. The details of the object modeling can be controlled upon vary requests. Meanwhile, these methods are easy to combine computer graphics with physical laws [Ang03, Wik07].

### 3.2 Procedural Methods

We use procedural modeling methods in our elastic object simulation system. One of the possible approaches to procedural modeling, a particle system, is designed to model elastic objects as described in this section. This particle system is also capable of describing the complex behaviors of elastic objects based on physical laws, such as solving differential equations of thousands of particles on those elastic objects. Another approach is language-based procedural method [Ang03], which generates complex objects with simple programs.

In order to model an elastic object, we need to study the following basic data structures, which are varied in one-dimensional, two-dimensional, and three-dimensional modeling methods.

###### Particles

are objects that have mass, position, velocity, and forces applied on them, but have no spatial extent. Moreover, the differential position and velocity change are two properties for these computation of each particle.

###### Springs

are massless with natural length not equal to zero. They are the linkage of particles. There must be at least one spring connects with two particles paired by modeling algorithm.

###### Faces

are the data type that consists of springs as the edges and particles as the vertices.

### 3.3 1d

The techniques used in an one-dimensional object are presented here, which are applied subsequently to models in two and three dimensions. An one-dimensional object with one end fixed as shown in Figure 3.1(a). The other end is interacted by users with mouse as in Figure 3.1(b). The interacted force is restricted to one dimension.

#### 3.3.1 Geometric Data Type

###### Particle

There are two mass particles and on a single spring .

###### Spring

In one-dimensional object, only one type of the spring, structural spring, is introduced. Structural spring, is used to model the object shape, connected by the two mass particles in this case. In Figure 3.1(a), the spring is at the initial state of equilibrium. The natural length of the spring is and the force is zero. In Figure 3.1(b), the spring is compressed by an external mouse force. The current spring length and the spring force restores the elastic object to its equilibrium position . When the spring is stretched out as shown in Figure 3.1(c), the current spring length and the force of the spring .

#### 3.3.2 Modeling Algorithm

• Step1: Create two particles and with positions and shown in Figure 3.2.

• Step2: Create a spring with these two particles as two ends and .

### 3.4 2d

In this section, we extend the one-dimensional elastic object to two dimensions. We create two separated elastic circles, inner circle and outer circle. Both of them consist of the same modeling structure as one-dimensional mass particles and springs. Then, the two concentric circles, inner and outer, are connected by various springs to become one two-layered elastic object. However, the distinct features in two-dimensional object have more types of springs presented and the air pressure inside the two-layer close shape is calculated. The spring surface prevents infinite expansion of the air; meanwhile the internal pressure avoids the surface collapsing.

#### 3.4.1 Geometric Data Type

Two-dimensional object is made of three types of primitives, mass particles, springs, and indexed triangular faces.

###### Particles

are defined based on their coordinates related to and axes. Consider a unit circle with twelve particles as an example shown in Figure 3.3. The spatial position for each particle is can be defined by the two equations:

where

is a small angle stepping along

###### Springs

In additional to the structural spring, which also exists in one-dimensional object, there are two other types of springs, radius spring and shear spring.

Structural springs: give the basic structure of inner circle and outer circle to prevent neighboring particles from getting too close to one another as shown in Figure 3.4(a). Linkage of each structural spring i is to connect with two particles and or and .

Radius springs: are the springs connected from particles on inner circle to the particle on the outer circle as part of the circle radius in order to prevent the bending of the surface as shown in Figure 3.4(b). Linkage of each radius spring is to connect with particle and the particle .

Shear springs: are springs connected from particles on inner circle to their diagonal neighbors on outer circle in order to avoid the surface fold over . Linkage of each left shear spring i is to connect with particle diagonally and the particle ; connect with particle diagonally and the particle and so on as shown in Figure 3.4(c). Linkage of each right shear spring i is to connect with particle diagonally and the particle ; connect with particle diagonally and the particle and so on as shown in Figure 3.4(d).

###### Faces

are the data structure for the only purpose of drawing and displaying a filled object to a two-dimensional object. The triangular facets can be drawn separately and visualized as a filled circle as shown in Figure 3.5.

#### 3.4.2 Modeling Algorithm

• Step 1: Define the number of particles as in our example. Then, the step size is degrees.

• Step 2: Define the group of particles’ position on inner circle and the ones on outer circle as shown in Figure 3.6. The first particle is at where , the second particle is at where … By multiplying the inner and outer coordinates with a different radius value, for example, , and to create two concentric circles.

• Step 3: Add the structural springs , , … to the inner circle according to the spring index of inner particles as shown in Figure 3.6. The same method is applied to outer structural springs on outer circle. The last spring, in our example, is composed of two particles and as two ends in order to make a close shape.

• Step 4: Loop through the number of structural springs to add the same number of radius springs according to the linkage of the inner and outer particles as shown in Figure 3.7.

• Step 5: Loop through the number of structural springs to add the same number of shear left springs and shear right springs according to the linkage of the inner and outer particles as shown in Figure 3.8.

### 3.5 3d

In this section, a more complicated three-dimensional elastic object is extended from the two-dimensional object. In the two-dimensional model, the structural springs’ index is the most important key data structure to link up all the particles and reference about the index of the particles. This spring linkage method will still work for the model based on the non-uniform sphere geometric modeling method. However, in the other geometric modeling method, the uniform sphere modeling, the faces’ index is the key data structure of the linkage to other data structure, such as particles and springs. The reason is because in the later geometric modeling method, each facet on the object is used for subdivision of other facets in each iteration. Compared with a two-dimensional object, the three-dimensional object consists the same types of primitives, such as particles, springs, and faces, but extended to axis.

#### 3.5.1 Non-Uniform Sphere

##### 3.5.1.1 Geometric Data Type

One of the simplest non-uniform modeling methods to generate an approximate facet sphere uses Polar to Cartesian Coordinates method. Consider the angle on -plane (around -axis), known as the Azimuthal Coordinate. The angle is from -axis, known as the Polar Coordinate. If we fix and draw curves as we change , we get circles of constant longitude; if we fix and vary , we obtain circles of constant latitude [Ang03].

###### Particles

The spherical coordinates for a particle can be defined by the three equations:

where

By stepping and in small angles and between their bounds as the number of slices and stacks, the particles are:

=

=

=

=

However, at the North and South Pole areas, we can only use triangles to present because all lines of latitude are converged.

The particle at the North Pole area can be presented as:

=

The particle at the South Pole area can be presented as:

=

###### Springs

There are also three types of springs in three-dimensional objects as we described in two-dimensional objects, such as structural, radius, and shear springs.

• Structural spring is still the basic data structure to form the shapes of inner and outer spheres. Four particles define four springs as the proper order. Taking the first four particles , , , and as an example, the first four springs’s are , , , shown in Figure 3.10. The structural springs on two poles are also defined by the particles on poles as proper order.

• Radius and shear springs, which connect inner and outer layers, follow the same methods as in two-dimensional object.

###### Faces

Any quadrilateral-facet on the body of sphere can be represented by four springs: , , , and . Any triangular-facet on the poles can be represented by three springs: , , and .

##### 3.5.1.2 Modeling Algorithm
• Step 1: Define the number of slices and stacks of a sphere, and in our example. Then, the step size is and .

• Step 2: Define the group of particles’ position on inner circle and the ones on outer circle shown in Figure 3.11. By multiplying the inner and outer coordinates with a different radius value, for example, , and to create two concentric spheres.

• Step 3: Add the structural springs , , … to the inner circle according to the spring index of inner particles as shown in Figure 3.11. The same method is applied to outer structural springs on outer circle.

• Step 4: Loop through the number of structural springs to add the same number of radius springs and shear springs according to the linkage of the inner and outer particles as described in two-dimensional object modeling method, shown in Figure 3.7 and Figure 3.8.

#### 3.5.2 Uniform Sphere

An important drawback of the non-uniform sphere model is that the faces vary in both shape (some are triangles and some are quadrilaterals) and size.

##### 3.5.2.1 Geometric Data Type

Surface refinement is a simple way for uniform modeling. It is started with a kernel polyhedron, which is a regular polyhedron with faces that are equilateral triangles. We have used an octahedron with bisecting each face at the same time recursively. This method is a powerful technique for generating approximations to curves and surfaces of a sphere to any desired level of accuracy.

###### Particles

The algorithm starts with a regular octahedron shown in Figure 3.12(a). The shape is composed of eight equilateral triangles, determined by six vertices, , , , , , and . The vertices of the kernel polyhedron are known to lie on the surface of a unit sphere of radius . We fix the two vertices and on axis and normalize the other five vertices by multiplying in order to make them lie on the unit sphere, centered at the origin. The six vertices after normalization are , , , , , and .

###### Faces

We talk about faces before talking about springs because the face is the key data structure for recursive subdivision and its index is referenced by spring index. The first eight triangular faces defined by the six particles are , , , , , , , .

###### Springs

Each face is composed of three springs. Therefore, the first twelve springs on the octahedron are , , , , , , , , , , , and .

###### Subdivision

We can subdivide a single triangular face of the kernel polyhedron by projecting the midpoints , , of its three edges onto the surface of the sphere as shown in Figure 3.13.

This face is split into four faces by bisecting each edge. The four new triangles are still in the same plane as the original triangle. We move the new vertices , , to the unit sphere by normalizing each new vertices. The number of particles increases by a factor of 2. The number of springs increases by a factor of 3. The number of facets increases by a factor of 4. We subdivide another 7 triangles with the same method. After subdividing the octahedron once, the number of particles are 12, the number of triangular faces is 32, and the number of springs is 36. We can repeat the subdivision process times to generate successively closer approximations to the sphere.

##### 3.5.2.2 Modeling Algorithm
• Step 1: Define a collection of particles to create a closed equilateral triangles shape of the elastic object. Define an octahedron object as the initial object with 6 particles, 8 triangular faces, and 12 structural springs as shown in Figure 3.14. Figure 3.14: Data Structure for Three-dimensional Uniform Object Representation without Subdivision Figure 3.15: Data Structure for Three-dimensional Uniform Object Representation with the Number of Subdivision n=1
• Step 2: Connect the particles with the structural springs according to the edge order of the octahedron to make an inner layer of the three-dimensional object. Each particle is separated equidistantly from its neighbors.

• Step 3: Check if there is need of subdivision to approach a more spherical object.

• Step 4: In the first subdivision, the object becomes 12 particles, 32 triangles, and 36 structural springs. The Figure 3.15 shows how new data can be inserted into a collection of particles, faces, and springs after subdivision.Use the first triangle as a concrete example. In the initial octahedron shape, find the midpoint of each edge on . Normalize the coordinates of these three new particles to make them lie on the sphere. Push these three new particles to the particle container. The first face of the initial octahedron has three pointers that point to particles , , and as shown in Figure 3.14. After subdividing this triangle, it becomes four smaller triangles connected by the bisectors , , and . The three new triangles are pushed onto container . The middle triangle replaces the original big triangle because the original triangle does not exist anymore. The pointers on each face point to the correspondent particles as indicated in Figure 3.15. New structural springs are added to Spring container correspondent to new faces only if there is no such spring has existed yet. The subdivision of the remaining faces follows the same method. More faces are approaching the object to a unit facet sphere.

• Step 5: Repeat Step 1 to 5 to create an outer layer with larger radius value than the inner layer.

• Step 6: Loop through the number of structural springs to add the same number of radius springs according to the linkage of the inner and outer particles.

• Step 7: Loop through the number of structural springs to add the same number of shear left springs and shear right springs according to the linkage of the inner and outer particles.

#### 3.5.3 Comparison of Non-uniform and Uniform Methods

The advantages of both methods are that they can be used to describe complex behaviors combined with physical laws, such as elasticity. Additionally, the level of detail (LoD) of the object can be adjusted depending on the proximity of the object on the display to the human’s eye.

The disadvantage of the non-uniform sphere modeling method is that the facets of the sphere do not have approximately equal size. The facets become smaller at the poles and bigger at the “equator”. Therefore, the springs are shorter at the poles and longer at the equator. The normal of each spring varies from equator and the ones on the poles. Consequently, this non-uniform modeling method increases errors in force computations for each particle.

The disadvantage of the uniform facet sphere generation algorithm is that it can not generate surfaces of arbitrary resolution. It can be shown that at all levels of recursion, particles at the kernel points are connected to four springs if the kernel object is an octahedron (as shown in Figure 3.12(b)). In other cases, all the particles at the kernel points are connected to five springs if the kernel object is an icosahedron (20 faces); all the particles at the kernel points are connected to three springs if the kernel object is a tetrahedron. All particles at recursively generated points are connected to six springs. This will result the irregular surface stiffness and might cause the non-spherical shape because the same pressure will displace the regions of a surface about the kernel points further than the rest of the surface.

To solve this problem, the sum of the spring forces accumulated at a particle can be normalized by multiplying a factor of , where is the number of springs connected to this particle. For example, if is a kernel point, which is connected to four springs and the sum of the spring forces is ; and if , which is the point generated from subdivision, connects to six springs and the sum of the six spring forces is . is multiplied by a factor of and is multiplied by a factor of .

Our simulation system ignores the described drawback resulting from the uniform sphere modeling method. We find a set of air pressure and spring stiffness parameter values at which the simulation is stable by trial and error. Thus, the difference of the forces for every particle either connected to four springs or six springs is not addressed in this work.

### 4.1 Gravity Force

Gravity force is a constant force at which the earth attracts objects based on their masses. In most cases, the particle system does not include gravitation, but, in our system, particle gravities represent object’s density. Users can set particle gravities to a non-zero value. g is a constant scalar of the gravitational field.

 Fg=mg (4.3)

### 4.2 Spring Hooke’s Force

Spring force is a linear force exerted by a compressed or stretched spring upon two connected particles. The particles which compress or stretch a spring are always acted upon by this spring force which restores them to their equilibrium positions. It is calculated as following according to Hooke’s law, which describes the opposing force exerted by a spring.

 Fh12=−(||r2−r1||−rl)ks (4.4)

where

is the first particle position,

is the second particle position,

is default length of the resting spring between the two particles,

is the stiffness of the spring,

when , the spring is resting,

when , the spring is extending,

when , the spring is contracting.

We have discussed the type of structural spring in one-dimensional object. In two and three dimensional object model, the same method applies on the other three types of spring, such as radius springs, shear left springs, and shear right springs with different spring stiffness and spring damping factor. So, the total Hooke’s spring force is:

### 4.3 Spring Damping Force

Spring damping force is also called viscous damping. It is opposite force of the Hook spring force in order to simulate the natural damping and resist the motion. It is also opposite to the velocity of the moving mass particle and is proportional to the velocity because the spring is not completely elastic and it absorbs some of the energy and tends to decrease the velocity of the mass particle attached to it. It is needed to simulate the natural damping due to the forces of friction. More importantly, it is useful to enhance numerical stability and is required for the model to be physically correct [BA97].

 Fd12=(v2−v1)⋅(r2−r1||r2−r1||)kd (4.6)

where

is the direction of the spring,

and is the velocity of the two masses,

is spring damping coefficient.

When the two endpoints moving away from each other, the force imparted from the damper will act against that motion; when the two endpoints moving toward each other, the damper will act against the squeeze motion. The damper will always acts against the motion. The total spring damping force is:

### 4.4 Drag Force

Drag force is the force when users interact with the elastic object through mouse. At the moment users click the mouse, the simulation system finds the nearest particle

to the current position of the mouse. If users drag this particle , the drag force contributes to force of this nearest particle. The forces applied on rest particles are effected by the new user applied force, which is passed through by springs.

We consider one end of the string connects to the mouse position and the other end of the string connects to the nearest particle on the object. This string is elastic, so it has all the spring’s properties, such as hook spring force and damping force. The drag force can be presented as following:

 Fa=−(||rm−ri||−rlm)ksm+(vm−vi)⋅(rm−ri||rm−ri||)kdm (4.8)

where

is the mouse position,

is the particle position nearest to mouse,

is default length of the resting mouse spring,

is the stiffness of the mouse spring,

is the velocity of the mouse represented as a mass

is the mass for the nearest particle,

is spring damping coefficient for the mouse spring.

is a momentary force for interacting with the elastic simulation system. This force is accumulated to the current forces already applied on this nearest particle.

In a one-dimensional object simulation system, the nearest particle to mouse is either or . In a two-dimensional and three-dimensional object system, the drag force is only applied on the outer layer of the double layered object when user interacts the object with mouse.

### 4.5 Air Pressure Force

In order to describe an elastic object more accurately, especially soft body of human beings and animals, the calculation only about the elastic force on the object’s surface is not enough. We add the flow pressure force inside of the elastic object to make the object wobbly looking when it is deformed.

The pressure force will be calculated for every spring, then update each particle’s direction. The pressure vector is always acting in a direction of normal vectors to the surface, so the shape will not deform completely. If pressure is simulated without also simulating the mass-spring system, the object will explode. (a) The External Gravity Force Is Applied Producing A Pressure Wave

In Figure 4.1(a), the object is deformed from bottom because of the gravity force. If there is no internal air pressure, the object will collapse unless the springs are hard enough to avoid the failure. With the very hard springs, it is difficult to simulate the reality of the elasticity. The real elastic object restores its shape as described in Figure 4.1(b). The simplified version of the Ideal Gas Law [Mat03], also known as Clausius Clapeyron Equation, is used to describe such effect:

 PV=NRT (4.9)

where

is the pressure value,

is the volume of the object,

is number of mols,

is the gas constant,

is the gas temperature. Therefore, the pressure force is:

 Fp=Pn (4.10)

where

is the pressure force vector,

is the normal vector to the springs on the object.

 P=NRTV (4.11)

#### 4.5.1 Volume

In order to find an estimate pressure inside of the object which will be applied to particles later, we need to calculate the volume of the object. The approximation of the volume is calculated with Gauss’ Theorem [Ker07]:

 V=∫∫∫vf(x,y,z)dxdydz⟺V=∫∫∫vf(x,y,z)dV (4.12)

where triple integrals of a function define a volume integral of an elastic sphere. Moreover, triple integrals can be transformed into surface double integrals over the boundary surface of a region if the three-dimensional object is closed shape by divergence theorem [Ker07]:

 V=∫∫∫VΔFdV⟺S=∫∫SFdS (4.13)

where

is a vector field,

is the object volume,

is the object surface.

Double integrals over a plane region may be transformed into line integrals by Green’s Theorem in the Plane:

 ∫∫SΔFdxdy⟺∫LFdL (4.14)

where is the object edge and is the length of the edge.

Therefore, a triple integrals function shown in Eq.4.12, which defines a volume integral of an elastic sphere, can be transformed to line integrals as shown in Eq.4.15.

 V≈∫LFdL (4.15)

We assume on the line, the vector field , the simplified integration of body volume is [Mat03, Ker07]:

 ∫LFdL=i=NUMS−1∑i=012(x1−x2)nxdL (4.16)

where

is the volume of the object,

is the absolute difference of the line (represents spring here) of the start and end particles at axis x,

is the normal vector to this line (spring) at axis x,

is the line’s (spring’s) length.

#### 4.5.2 Normals

Normals are unit vectors perpendicular to specified data structure, such as particle (vertex psuedo-normals), spring (line), and face (polygonal facets) on the object.

• Particle normal, or vertex psuedo-normals, does not exist for vertices; however, it can be considered as the average of the normals of the subtended neighbor particles. To calculate the particle normal is to sum up the normals for each face adjoining this particle, and then to normalize the sum.

• Spring (line) normal in two dimension is based on the two particles connected on the spring. It is perpendicular to the spring itself.

• Face (plane) normal in three dimension is determined by right-hand rule, which is perpendicular to its surface based on the any pair of springs on the surface. The normal for a triangle surface composed with three particles is computed as the vector cross product of the springs and .

The usage of normal calculation method in our elastic object simulation system is for analysis of the direction of the pressure force inside of the object. Only spring normal is calculated here because all the internal and external forces will apply on each spring, and the spring will define the particles’ force which connect onto it.

##### 4.5.2.1 2D Normals

For the single spring , the Cartesian coordinates for particle is ; the Cartesian coordinates for particle is . The normal to this spring is the spring rotate at axis z according to the space position. So, we can get the components of the normal in x axis and y axis as following

 [x′y′]=[cos900−sin900sin900cos900][x2−x1y2−y1]
 [x′y′]=[−(y2−y1)x2−x1]
##### 4.5.2.2 3D Normals

The calculation of the 3D normals of springs is important because it will define the direction of the internal air pressure force, either compress the elastic object or expand its volume. In theory, in three-dimensional simulation, the normal of a spring in the space position is represented as an average of the normals of faces connected to it. However, in our elastic object simulation system, we use a simplified estimated normal based on the normal algorithm of the two-dimensional calculation discussed above. Instead of rotating a line at z axis to get its normal vector in two-dimension, the estimated algorithm is rotating a line at z axis, y axis, and x axis to get its normal in three-dimension.

 ⎡⎢ ⎢ ⎢⎣x′y′z′1⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢⎣cos900−sin90000sin900cos9000000100001⎤⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢ ⎢⎣cos9000−sin90000100sin9000cos90000001⎤⎥ ⎥ ⎥ ⎥⎦
 ⎡⎢ ⎢ ⎢ ⎢⎣10000cos900sin90000−sin900cos90000001⎤⎥ ⎥ ⎥ ⎥⎦⎡⎢ ⎢ ⎢⎣x2−x1y2−y1z2−z11⎤⎥ ⎥ ⎥⎦

Therefore,

 ⎡⎢ ⎢ ⎢⎣x′y′z′1⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢⎣z2−z1y2−y1−(x2−x1)1⎤⎥ ⎥ ⎥ ⎥⎦

We use a vector as an example to prove this algorithm, the normal vector for this vector is:

 ⎡⎢ ⎢ ⎢⎣x′y′z′1⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢⎣0−00−0−(1−0)1⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢⎣00−11⎤⎥ ⎥ ⎥⎦

This result is reasonably correct and believable despite of the fact that if vector lies in the -plane or lies in the -plane.

However, it shows that this estimation algorithm has the limitation for some cases, for example vector , the normal vector is:

 ⎡⎢ ⎢ ⎢⎣x′y′z′1⎤⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢ ⎢⎣0−01−0−(0−0)1⎤⎥ ⎥ ⎥ ⎥⎦=⎡⎢ ⎢ ⎢⎣0101⎤⎥ ⎥ ⎥⎦

This result shows the normal vector is the vector itself, which is obviously wrong. However, with this estimated algorithm, the simulation result appears enough realistic; moreover, it requires less computational effort111Our estimation only takes 3 additions vs. 12 multiplications and 6 additions for two cross products and three more additions and divisions for the averaging..

### 4.6 Collision Force

If an object continues traveling under a force without colliding with other objects, it will be very difficult to describe objects’ motion and elastic response in reality. Collision force is the force to make object bounce away from the fixed interacting plane when elastic object collision happens.There are two steps to describe the collision effects: detection and reaction. Detect the elastic object if particles hit anything; adjust their position by computing the impulse.

#### 4.6.1 Collision Detection

Collision Detection is a geometric problem of determining if a moving object intersected with other objects at some point between an initial and final configuration. In our elastic object simulation system, we are concerned with the problem of determining if any of particles collide with any of solid planes.

###### Perfect Elastic Collision

We take one particle collides with a plane shown in Figure 4.2 as an example. We can detect this collision by inserting the particle position into the plane equation:

 P(x,y,z)=ax+by+cz+d (4.17)

If , the particle is within the plane. If , the particle collides with the plane. If , the particle penetrates the plane. At each time step, looping through all the particles on the object, each particle is checked if it is outside of the interacting plane.

When the particle collides with the plane, if there is a perfect elastic collision as in Figure 4.2, the particle does not lose its energy, so its speed does not change. However, its direction after the collision is in the direction of a perfect reflection.

 Fc=2((P−P′)⋅n)n−(P−P′) (4.18)

where

is the the direction of a perfect reflection

is the normal at the point of collision and the previous position of particle

is the vector from the particle to the surface.

###### Damped Elastic Collision

If there is a damped elastic collision, the particle cannot penetrate the surface, and it cannot bounce from the surface because of the force being applied to it, then we need to apply the damped elastic collision method. The particle loses some of its energy when it collides with another object. The coefficient of restitution of a particle is the friction of the normal velocity retained after the collision. Therefore, the angle of reflection is computed as for the inelastic collision, and the normal component of the velocity is reduced by the coefficient of restitution.

#### 4.6.2 Collision Response

Collision Response is a physics problem of determining the forces of the collision. In elastic collision, elastic object should bounce away from the colliding plane and some energy is lost in the collision response as described in the penalty method.

 Fc=−eFc (4.19)

where is elasticity of the collision and . At , the particle does not bounce at all; , the particle bounces with no friction.

In an one-dimensional object, the boundaries are the walls and floors. In a two-dimensional and three-dimensional object, the particles on the outer layer still follow the same method and same pre-defined boundary as the one-dimensional object. However, for the particles on the inner layer, the boundary is constrained to the outer layer instead of the wall and floor.

### 4.7 Force Accumulation Algorithm

The following algorithm describes how different forces are accumulated and applied to an elastic object. For a one-dimensional object, some steps will be skipped, for example, there are no other types of spring computations except structural springs because other types of springs only apply on two-layer 2D or 3D objects. Moreover, there are no pressure force accumulation and volume computation because these steps are only available for closed shape objects.

• Step 1: Loop through the number of particles to assign particles with mass value and compute gravity force . Gravity force, which is independently on each particle, either depends on a constant force, or one or more of particle position, particle velocity, and time [Wit97]. If the object is one-dimensional, the mass of each particle can be different. If the object is two or three dimensional, the mass of the particles on inner or outer layer can also be set differently.

• Step 2: Loop through the number of the structural springs to accumulate the structural spring force.

• Step 3: Loop through the number of the radius springs to accumulate the radius spring force.

• Step 4: Loop through the number of the shear springs to accumulate the shear spring force.

• Step 5: Initialize density as gas, liquid, or rubber inside of the body and introduce some simple physics to describe it. In the current system, only air pressure material is considered and only pressure equation will be used for this extra force computation.

• Step 6: Calculate volume of the inner layer and outer layer of the elastic object.

• Step 7: Calculate the normals of springs on each triangular face to define the pressure force direction.

• Step 8: Calculate the force from the internal air pressure by multiplying the force value by normal vector of the spring.

• Step 9: Accumulate pressure force to each particle.

• Step 10: If users apply the drag force, compute the user applied force and accumulate this force to the dragged particle.

• Step 11: Integrate the object’s momentum motion by calculating the derived velocity and its new position for each particle. This step will be explained in next chapter.

• Step 12: Resolve collision detection and response and define the updated position.

### 5.1 Differential Equations

Differential equations describe a relation between a function and one or more of its derivatives. The order of the equation is the order of the highest derivative it contains. The elastic object simulation system is associated with initial value problems because it always seeks the particles’ velocity and position at next time step from their initial state at time . We will concentrate on ODE (ordinary differential equation), where all derivatives are with respect to single independent variable, often representing time, such as position and velocity, during the derivate of the state at discrete time steps [Ang03].

 y′=A(y,t) (5.1)

where

is a function of and ,

is a vector, which is the state of the system,

is a vector, which is ’s time derivative.

Suppose that we integrate the Eq.5.1 over a short time

 y(t+h)−y(t)=∫t+htA(y,t)dt (5.2)

where

is the small stepsize of time,

is the initial state at the start point ,

is the value we need to find over time thereafter.

Thus

 y(t+h)≈y(t)+hA(y(t),t) (5.3)

#### 5.1.1 Explicit Euler Integrator

The simplest ODE integration method is Explicit Euler Integration method or Forward Euler method. It evaluates the forces at time , compute derivatives at the state of by multiplying the interval , and add it to the current state . Consider a Taylor series expansion as in Eq.5.4:

 y(t+h)=y(t)+hy′(t)+h22!y′′(t)+h33!y′′′(t)+⋯+hnn!(∂ny∂tn)+⋯ (5.4)

Euler method retains only first derivative:

 y(t+h)=y(t)+hy′(t)+O(h2) (5.5)

We split the series into elements, which we will later use in a re-usable manner throughout integrator framework, where

which represents the first term in Eq.5.5, is the initial state

 k0=y(t) (5.6)

which represents the second term in Eq.5.5, is the function to find the simplest estimation, the Euler slope of the interval.

 k1=y′(t)=A(y(t),t) (5.7)

Thus

 y(t+h)=k0+hk1 (5.8)

We can apply this method iteratively to compute further values at state , ,…. [BD03] This method is easy to implement; however, it is a low accuracy prototype ODE. In Figure 5.2, we can see Euler method only calculates the derivative, also called slope, at the beginning of the interval and adds it to the value at the initial state; therefore, it is asymmetric and not stable. .

###### Pseudocode for Euler Method

Line 1: define A(y(t), t) Line 2: initial values y0 and t0 Line 3: stepsize h and number of steps n Line 4: for i from 1 to n do Line 5: k1 = A(y(t), t) Line 6: y = y + hk1 Line 7: t = t + h

#### 5.1.2 Midpoint Integrator

Compared to the Euler method, the one-sided estimate algorithm, midpoint integrator is a symmetric estimate method with a higher per-step accuracy. It computes the derivative at the center of the interval first, then computes the end of the interval.

The midpoint integrator, just like others, is based on the Taylor’s series. It retains only first three derivative term:

 y(t+h)=y(t)+hy′(t)+h22!y′(t)+O(h3) (5.9)

We split the series into elements again for explanation of the method, where

, which represents the first term in Eq.5.9, is the initial state at time .

 k0=y(t) (5.10)

which represents the second term in Eq.5.9, is the function to find the the simplest Euler slope of the interval at time .

 k1=y′(t)=A(y(t),t) (5.11)

is the function to find the the simplest Euler slope of the interval at time .

 k2=y′(t+h)=A(y(t+h),t+h) (5.12)

Since the unknown appears on the right side of Eq.5.13, in as one of the arguments of function , we can use the value obtained using the Euler method in Eq.5.5.

 A(y(t+h),t+h)≈A(y(t)+h(y(t),t),t+h)=A(y(t)+hk1,t+h) (5.13)

The midpoint integration technique obtains a more accurate estimate of the slope than Euler’s technique. The following equation computes the integrand at the middle of the interval of and shown in Figure 5.3. Thus,

 y(t+h)=k0+hk1+k22 (5.14)

Compared to Euler Method, Midpoint Method, also called the Runge-Kutta method of order 2, goes from to , we must evaluate function twice. By using Taylor’s theorem to evaluate the per-step error, we would find that it is now . Therefore, this method is more stable than Euler Method with same step size.

###### Pseudocode for Midpoint Method

Line 1: define A(y(t), t) Line 2: initial values y0 and t0 Line 3: stepsize h and number of steps n Line 4: for i from 1 to n do Line 5: k1 = A(y(t), t) Line 6: k2 = A(y(t+h), t+h)= A(y+hk1, t+h) Line 7: y = y + h/2(k1+k2) Line 8: t = t + h

#### 5.1.3 Runge Kutta Fourth Order Integrator

Runge Kutta Fourth Order integrator evaluates the derivative four times. It is the most accurate integrator that we describe compared to Euler and Midpoint.

The Runge Kutta Fourth integrator, is also based on the Taylor’s series. It retains only first five derivative term with a local truncation error :

 y(t+h)=y(t)+hy′(t)+h22!y′′(t)+h33!y′′′(t)+h44!y′′′′(t)+O(h5) (5.15)
 k0=y(t) (5.16)
 k1=y′(t)=A(y(t),t) (5.17)
 k2=A(y(t)+hk12,t+h2) (5.18)
 k3=A(y(t)+hk22,t+h2) (5.19)
 k4=A(y(t)+hk3,t+h) (5.20)
 y(t+h)=k0+16h(k1+2k2+2k3+k4) (5.21)

where

is the initial state

is the slope at the left end of interval,

is the slope at the middle point using the Euler formula to go from to ,

is the second approximation to the slope at the midpoint,

is the slope at using the Euler formula and the slope to go from to .

###### Pseudocode for Runge Kutta Fourth Order Method

Line 1: define A(y(t), t) Line 2: initial values y0 and t0 Line 3: stepsize h and number of steps n Line 4: for i from 1 to n do Line 5: k1 = A(y(t), t) Line 6: k2 = A(y+h/2(k1), t+h/2) Line 7: k3 = A(y+h/2(k2), t+h/2) Line 8: k4 = A(y+hk3, t+h) Line 9: y = y + h/6(k1+2*k2+2*k3+k4) Line 10: t = t + h

### 5.2 Newton’s Laws

After the force accumulation on the object, it is important to find the acceleration a in order to define the motion of objects in their next time step. The physical law that governs the motion of objects is the Newton’s Second law. It states that the force is proportional to the time rate of change of its linear momentum. Momentum is the product of mass and velocity v.

 F≈mΔvΔt (5.22)
###### Velocity

v is the integral of acceleration a with respect to the time . Therefore, integrating the acceleration gives us the new velocity v.

###### Position

r is the integral of velocity v with respect to the time . Therefore, integrating the velocity gives us the new position r.

 r=∫vdt (5.24)

Let’s take one particle on the object as an example and understand how the different integrators work.

#### 5.2.1 Newton’s Laws in Euler Integrator

Based on the Euler Integrator method shown in Eq.5.5, the new velocity and position of a particle can be integrated follows.

###### Velocity

can be represented as the following equation:

 v(t+h)≈v(t)+hv′(t) (5.25)

represents the first term in Eq.5.25, which is the initial velocity at time

 vk0=v(t) (5.26)

represents the second term in Eq.5.25, which is the function to compute the derivative velocity in the period

 vk1=hv′(t)=a(t)h (5.27)
###### Position

can be represented as the following equation

 r(t+h)≈r(t)+hr′(t) (5.28)

is the initial position at time

 rk0=r(t) (5.29)

is the function to find the travel position in the period

 rk1=hr′(t)=v(t)h (5.30)

#### 5.2.2 Newton’s Laws in Midpoint Integrator

We apply the midpoint algorithm theory on the Newton’s law in order to achieve higher accuracy in the the relationship between the velocity and the position according the Eq.5.9.

###### Velocity

can be represented as the following equation

 v(t+h)≈v(t)+hv′(t)+h22!v′′(t) (5.31)

is the initial velocity at state

 vk0=v(t) (5.32)

is the function to compute the derivative velocity in the period

 vk1=v′(t)=a(t)h (5.33)

is the function to compute the derivative velocity in the period

 vk2=v′(t+h)=v(t)+a(t)h (5.34)

Therefore, the new velocity of a particle is

 v(t+h)=vk0+vk1+vk22 (5.35)
###### Position

can be represented as the following equation

 r(t+h)≈r(t)+hr′(t)+h22!r′′(t) (5.36)

is the initial position at state

 rk0=r(t) (5.37)

is the function to find the travel position in the period

 rk1=r′(t)=v(t)h (5.38)

is the function to find the travel position in the period

 rk2=r′(t+h)=r(t)+v(t)h (5.39)

Therefore, the new position of a particle is

 r(t+h)=rk0+rk1+rk22 (5.40)

#### 5.2.3 Newton’s Laws in the Runge Kutta Fourth Order Integrator

Based on the Runge Kutta Fourth Order method we have shown in Eq.5.9, the new velocity and position of a particle can be integrated as following.

###### Velocity

can be represented as the following equation

 v(t+h)≈v(t)+hv′(t)+h22!v′′(t)+h33!v′′′(t)+h44!v′′′′(t) (5.41)

is the initial velocity at time

 vk0=v(t) (5.42)

is the function to compute the derivative velocity in the period

 vk1=a(t)h (5.43)

is the function to compute the derivative velocity of the Euler integration in the period based on the previous step

 vk2=vk0+vk12 (5.44)

is the function to compute the derivative velocity of the second approximation based on the in the period

 vk3=vk0+vk22 (5.45)

is the function to compute the final resulting velocity change of from

 vk4=vk0+vk3 (5.46)

Therefore, the new velocity of the particle is

 v(t+h)=vk0+16h(vk1+2vk2+2vk3+vk4) (5.47)

If we integrate the velocity vector over time, it gives us how the position vector changed over this time.

###### Position

can be represented as the following equation

 r(t+h)≈r(t)+hr′(t)+h22!r′′(t)+h33!r′′′(t)+h44!r′′′′(t) (5.48)

is the initial position at time

 rk0=r(t) (5.49)

is the function to find the travel position in the period

 rk1=v(t)h (5.50)

is the function to find the travel position of the Euler integration in the period based on the previous step

 rk2=rk0+rk12=r(t)+v(t)h2 (5.51)

is the function to find the travel position of the second approximation based on the in the period

 rk3=rk0+rk22 (5.52)

is the function to find the travel position change of from

 rk4=rk0+rk3 (5.53)

Therefore, the new position of the particle is

 r(t+h)=rk0+16h(rk1+2rk2+2rk3+rk4) (5.54)

### 5.3 Comparison of Three Integrators

#### 5.3.1 Efficiency

For a given step size, Euler is more efficient because it requires only one derivative evaluation per step. Mid Point requires about twice as much computation than the Euler integrator because Mid Point uses two steps to calculate velocity and position at the next time. Runge Kutta Fourth Order requires about four times as much computation as Euler integrator because it use four steps to calculate the velocity and position at the next time step [BD03]. For some configuration, if speed is the priority, Euler integration is convenient to use, but at the expense of accuracy and stability.