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.
What form should deriv take when returned by evalODEFuntion? We know its an anonymous function, but the sizes in StepScene are not matching up.
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.
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!
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.
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?
@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).
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.
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.
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.
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?
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.
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.
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?
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.
@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.
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?
@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).
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.
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))]);
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
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?
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?)
@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?
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.
@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.
@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.
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?
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?
@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.
@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.
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?
@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.
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.
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.
@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.
@Daniel: Yes, you’re right, the size of the matrix is (# of springs) x 5.
Does matlab have a function to return the slope of the linear regression line of two data sets?
@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.
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?
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,
@mt438: At timesteps of those sizes forward euler may indeed explode for the bike scene.
@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.