Lab 8

From CSE100B Wiki

Differential Equations

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.

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

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.

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.

Student Questions and Answers

Questions? Please post them on the Google Group http://groups.google.com/group/wustl-matlab.