## Project 3

Project 3 is now posted. In this assignment you will build a simple mass-spring simulator that allows the user to load scene descriptions from a file, select different integrators, and interact with the simulation using both “probe” (a circular wall centered at the current mouse point) and “pull” (a spring between a given particle and the current mouse point) tools. It will also models gravity, viscous damping, spring, and wall collision forces.

### 43 Responses to “ Project 3 ”

• fml666

What form should deriv take when returned by evalODEFuntion? We know its an anonymous function, but the sizes in StepScene are not matching up.

• mev25

Deriv must match the size of the input parameter state. This is intuitively true as it should store the derivative of every component in state.

• student

Could you please explain the input parameters of odeForwardEuler? I imagine that tspan is the span of time? (a vector of length 2 containing the initial time and ending time?). Is u0 just a number representing the initial value? what about options? Which one tells me or passes in the step size?

Thanks!

• mev25

You can find the relevant information in section 3.7 but the quick answer to your question is that u0 represents the initial state of the system of ODEs (if you only had one differential equation then it would be a value but this is not the case). The step size is given by options.InitialStep.

• kslice

What is the correct format when using odefun? Can I do something along the lines of odefun(U,options), or do I have to convert something here?

• Manolis Savva

@kslice: odefun is a function handle with the signature odefun(t,u). In other words, it takes the current time and current state vector as arguments. You can create this function handle by using an anonymous function declaration as described in section 3.7 of the handout: @(t, u) evalODEFunction(u, scene).

• mev25

You can also find the relevant info on odefun in section 3.7 but basically, odefun takes a scalar time and a state vector u (because it represents G(t,u) in the writeup). Note that G only returns the derivative of the state at a particular time and thus does not depend on the timestep used by the integrator nor the precision, so there’s no reason to pass options as an argument.

• fk

It says in the writeup to treat walls as springs (when the walls are enabled). What is the “desired distance” r for these springs? If I set it to zero, the objects gets accelerated A LOT when they cross over walls.

Also, these ‘wall springs’ are only one-directional right? (Ie, they aren’t always pulling objects towards them?)

Thanks.

• mev25

The desired distance is in fact zero, otherwise the system would strive to keep objects inside the wall. You should only apply wall forces if and only if the object is outside bounds. If you are noticing objects accelerate significantly when they strike a wall then a small bug is likely.

• manischewitz

I still don’t understand how to access the variable ‘scene’ from odefun ( @(t, u) evalODEFunction(u, scene) ). I’ve tried the dot operator and other tricks but to no avail. How do I do this in Matlab?

• mev25

I think you are confusing odefun and evalODEFunction. evalODEFunction is a function that evaluates the derivative of a particular state in our springies environment and, as such, needs the scene information. If you’re writing a black box integrator, such as the ones that already exist in MATLAB and the ones we’ve asked you to write, you want it to work as general as possible which means that all you need is a black box function (known only to the integrator as odefun) that returns the derivative of the state. We express odefun in terms of evalOODEFunction which we shortcut by using anonymous functions (we have already done this for you).

TLDR: In the context of the integrators you have to write, forget about springies. Write general purpose integrators.

• sac76

I don’t understand exactly what we are supposed to do for the “questions” part (where we have to create the graph of the errors). The MAT files that we get from recording contain point data and time data, and you say we are supposed to use the length of the pendulum as an error metric. So for a given point in the plot (e.g. Forward Euler, timestep 0.001, precision 10^-5), are we supposed to look at PointData(:,:,1) corresponding to that session and use the distance between the two points as the value for the y-axis and Time(1) as the value for the x-axis?

Thanks.

• sac76

Sorry, a minor correction to my post above. Do you mean that we subtract the length of the pendulum in the ODE 45 data from the length of the pendulum in the Forward Euler data and use the absolute value of that as the value for the y-axis?

• sac76

Actually, even with the above correction my post still doesn’t make sense. Back to my original question, I don’t understand exactly what the graph is supposed to be of.

• Manolis Savva

@sac76: When calculating the absolute error in the length of the pendulum you need to compare the length you obtain with forward euler or midpoint against a more accurate method such as Matlab’s ode45. So you essentially only need to look at this difference for the first frame and then vary the step size to see how it affects the error. Your graph will be one of absolute error in pendulum length (with ode45 as a base-line) against the step size (i.e. h). So you were right in your second post.

• aip23

The handout refers to gravity in different parts as a force and as an acceleration. Just to be clear: is scene.Gravity.Direction the force (i.e. it needs to be divided by mass to get acceleration), or the acceleration (i.e. it is independent of mass) on a particle?

• Manolis Savva

@aip23: It is the acceleration and should be independent of mass (indicated by the fact that a single vector is defined for the entire scene — typically 10 m/s^2 in the negative y direction).

• Manolis Savva

In response to some questions about the expected performance of evalODEFunction: Performance will vary significantly with the particular hardware configuration. However, as a general example to illustrate what you should expect, running the bike scene with ode113 and the default timestep and precision on reasonable hardware you should get at least one frame per second. Simpler scenes should be much faster, generally running at interactive rates.

• mev25

If you would like to add a frames per second counter to your springies, add the following line (as a single line, since the blog wraps it) to the end of the mspringiesPainter.m file:

set(handles.figure1, ‘NAME’, ['FPS: ',num2str(floor(1/handles.timer.AveragePeriod))]);

• mev25

If you have a particularly fast implementation you can also reduce the max 30fps restriction that is currently in place by reducing the period on line 75 of mspringies.m (to like 0.005). This will make your animation smoother. :D

• student

I am implementing acceleration as a result of wall and I have trouble with the formula Ks(||Xj-Xi||-r)*(Xj-Xi)/||Xj-Xi||. If r=0, the ||Xj-Xi|| cancels out and we are left with a scalar mulfiplied by a vector, which would yield a vector? Are we then splitting it into the x-component and the y-component, and adding that to the Xacc and Yacc respectively?

• student

Hi,
Please ignore the post above I think I figured it out.
I do have a new question though… I am trying to calculate the acceleration due to wall damping… it says to take the difference in velocities and taking the dot produce with the direction of the spring… one of my point should be a point on the wall, and I don’t know what I should put for the velocity. Is the direction of the spring the direction the spring is pulling? (so -y for ceiling?)

• Manolis Savva

@student: I’m not sure if I understand your question correctly but if you are asking what velocity should be used in the expression for the virtual spring force that a wall applies to a point moving outwards then the answer is obviously the velocity of the particle itself, since the virtual point on the wall is stationary. The direction of the spring is something you define by the choice of which side you consider to be primary (if you look at the expression for the spring force, the normalized vector d points from x_i to x_j so its direction depends on which point you call x_i and which x_j).

When testing using Bike, the front wheel grows huge and the back wheel ends up disconnecting and becoming more of a line than a circle (this is with the stepsize and precision turned to the lowest settings). We were wondering what the ideal behavior is supposed to be?

• student

When it comes to wall forces, there is one with coefficient wallstiffness and one with coefficient WallDamping. I think I know how to do the first one but I am unsure of how to implement the second one… The handout mentions that damping of springs has a formula fi= kd((xj* – xi*)*d)d. The * next to xj and xi represent dots on them (velocity). So xj*-xi* is the difference between two velocities. One is that of the particle, is the other 0 because the point on the wall is stationary?

Thanks.

• Manolis Savva

@adm43: The bike should roll to the left side of the screen with its wheels rotating around their central axis. In other words you shouldn’t observe any “explosive” or “implosive” behavior like what you describe.

• Manolis Savva

@student: You are right that the second part is the damping force which will be the same as the damping force for a regular spring, except that the point on the wall has 0 velocity. So you can just evaluate the expression you quoted.

• Michael

When writing our odeMidpoint and odeEuler function, how do we use odefun? For example, in the Euler method, we are supposed to multiply the result of the derivative by the timestep. To get the result of the derivative, we need to evaluate our odefun. How do we get the result out of it?

• student

When calculating forces due to probe, is the rest length 0 as in the case with walls? I imagine that radius of probe is given so that we could find the closes distance between a point and the probe by finding the closest point on the circumference of the probe 2pi*r? Is that correct?

• Manolis Savva

@Michael: Since odefun is a function handle defined to conform to G(t,u) as described in the handout you would need to use it like so: odefun(t,u) where t is the current time and u is the current state column vector. Refer to section 3.7 of the handout for more details.

• Manolis Savva

@student: Yes, the probe behaves just like the regular walls except it has a circular shape. So as soon as you determine the closest point on the circle to a given particle you apply exactly the same approach as for the wall collisions.

• mt438

To use midpoint euler, we would have a term like f(t+h/2,u+k1/2). Because the force is independent of time, so I think odefun(t, u+k1/2) and odefun(t+h/2, u+k1/2) will return the same result. Is it corret?

• Manolis Savva

@mt438: Yes, you’re right. All the forces we use do not have an explicit dependence on time. However, in general, this may not be the case which is one reason for odefun to take time as its first argument.

• fk

there should only be one type of damping on a particle at any point in time, right? Because if a particle is in a wall (or in the probe) it can’t also be in the atmosphere of the rest of the setting.

Thanks.

• Daniel Brooks djb332

I’m pretty sure there’s a typo in the handout. It says that the size of the spring matrix is n by 5 matrix where n is the number of particles. This is incorrect, as spring is a (# of springs) x 5 matrix.

• Manolis Savva

@fk: There may be more than one type of damping on a particle since we have viscous damping for all particles with a non-zero velocity and also spring force damping for particles that are attached to a spring. You don’t have to worry about “disabling” viscous damping outside the screen region.

• Manolis Savva

@Daniel: Yes, you’re right, the size of the matrix is (# of springs) x 5.

• mt438

Does matlab have a function to return the slope of the linear regression line of two data sets?

• Manolis Savva

@mt438: You can use matlab’s curve fitting functionality available from Tools menu, Basic Fitting at the top of a Figure window. This only works on one data set at a time. It should also be relatively easy to manually measure the gradient by taking the typical “dy/dx” approach.

• mt438

I observe some explosive behavior when I try to run scenes such as the bike with forward euler method and timestep = 0.25. Other methods work well with the same setting. Is this expected?

• student

Hi,
What are you supposed to submit for this project? I zipped the m files I actually coded but the handout says to submit “all MATLAB code you used.” Does this include other files included in the framework?

Thanks,

• Manolis Savva

@mt438: At timesteps of those sizes forward euler may indeed explode for the bike scene.

• Manolis Savva

@student: You can zip only the .m files you actually changed if you want. Be sure to include some documentation answering the required questions as well.