From CSE100B Wiki
Systems have state variables such as the position of the puck on the hockey table (Lab 6 and 7). Many times, the state at time t can be determined from the state at time (t - 1), the state at time (t - 2), etc. For example, given a constant velocity, the position of the hockey puck at time t was given by: position(t) = position(t - dt) + dt * velocity. Given an initial position and velocity, this equation allows us to calculate the position of a puck for times in the future (and times in the past). This equation is an example of a differential equation. Differential equations (DE) relate the state variables to changes in the state variables. Modeling a system as a set of differential equations is a powerful technique and MATLAB contains many solvers for simulating these systems. We could model the hockey table by a set of differential equations that relate position, velocity, and acceleration due to the walls, repelling forces, etc., but it would be messy.
To gain familiarity with solving DEs with MATLAB, we will use MATLAB to simulate a simple spring-mass system. We'll start by simulating it like the hockey table.
Q1: Read the following script which contains a simulation of a mass-spring system.
m = 3; % Mass k = 8; % Spring constant y = -1; % y position vy = 0; % y velocity ay = 0; % y acceleration g = 9.80665; % Acceleration due to gravity dt = .1; % Time step y_0 = 0; % Resting position of mass with no gravity. % Draw the initial state. line ([-2 2], [5 5]); hold on axis ([-2 2 -10 6]); mass_handle = plot (0, y, 'og', 'MarkerSize', 10); spring_handle = line ([0 0], [5 y], 'LineStyle', '--'); for step = 1:1000 % Calculate the force from the spring. F_spring = -k * (y - y_0); % Calculate the force from gravity. F_gravity = -m * g; % Calculate the acceleration due to the forces. ay = (F_spring + F_gravity) / m; % Update the velocity. vy = vy + dt * ay; % Update the position. y = y + dt * vy; % Draw current state. set (mass_handle, 'YData', y); set (spring_handle, 'YData', [5 y]); pause (.01); end
Run the script and observe how the mass moves. Estimate the maximum and minimum heights of the mass during the simulation.
To use MATLABs ODE solvers, we must convert the system to differential equations. Specifically, we need to create a function that takes the current simulation time and state and returns the updated state. Our function will take a column vector where the first element is the position of the mass and the second element is the velocity. Our function will return a column vector where the first element is the velocity and the second element is the acceleration. This function is shown below.
function change_in_states = mass_spring_update(time, states) y = states(1); vy = states(2); m = 3; % Mass k = 8; % Spring constant g = 9.80665; % Acceleration due to gravity y_0 = 0; % Resting position of mass with no gravity. % Calculate the force from the spring. F_spring = -k * (y - y_0); % Calculate the force from gravity. F_gravity = -m * g; % Calculate the acceleration due to the forces. ay = (F_spring + F_gravity) / m; change_in_states = [vy; ay];
Notice that most of the lines are take directly from the script above.
To perform the simulation we use the following script:
y = -1; % Initial y position vy = 0; % Initial y velocity states = [y; vy]; % Pack into a column vector time = [0 10]; % Specify the start time and end time [time_vector pos_vel_vector] = ode45(@mass_spring_update, time, states); % Simulate !! figure; plot(time_vector, pos_vel_vector(:, 1)); % Plot the position of the mass vs. time
The @mass_spring_update creates a reference to a function. This prevents MATLAB from evaluating it before the call to ode45 which would only pass its output to ode45. Instead, it allows MATLAB to "remember" the function and call it later.
Q2: Copy-and-paste the function. Save it as mass_spring_update.m. Remember, MATLAB calls the function by the name of the file. Copy-and-paste the script. Save it as lab8_1.m because you will modify it in the future. Run lab8_1. What is the maximum and minimum position of the mass? What is the maximum and minimum velocity of the mass? How many time steps did MATLAB use to solve the differential equation? (Don't estimate. Use MATLAB to calculate these answers.)
lab8_1.m: Add damping to the mass-spring system. The equation for the force of damping on the mass is F = -c * v where c is the damping constant and v is the velocity. Pick a spring constant and simulate the system long enough to see the mass stop moving.
Often, we need to simulate systems using different values of a parameter such as the damping constant. Currently, this constant is defined inside our function. If we want to change it from outside of the function, we must declare it as a global variable.
lab8_2.m: Plot the position of the mass for a range of different damping constants.
Demonstrate Q1-Q2 and lab8_1, lab8_2.m to the TA or instructor.
Student Questions and Answers
Questions? Please post them on the Google Group http://groups.google.com/group/wustl-matlab.