Deliverables
You will upload to Canvas a single .zip file, titled LastNameFirstInit_Lab8.zip, containing two separate files: 1) LastNameFirstInit_Lab8.m will contain your ABM3 code in response to Task 1, and 2) LastNameFirstInit_Lab8.pdf will contain your lab report in response to Task 2.

Terminology
State Variables
In a dynamical system (such as an Initial Value Problem), the dependent variables are also called the state variables. e.g., in Problem 5 of HW #9, x and v (or ẋ) are the state variables. These variables are often put into an array, which is called the state array. While the independent variable (usually time, for us) is not technically considered a state variable, we are primarily interested in the system state as a function of time. In our case, we will insert the independent variable as the first value in the state array and call it the augmented state array.

Mode Shapes and Natural Frequencies
Consider the vibration system below. As a two degree of freedom system, it has two natural frequencies and two mode shapes. If the two masses are both displaced to the left (by the correct distances) and released, they will each undergo perfect sinusoidal motion. Both will oscillate at the same frequency – called the first natural frequency. The first mode shape is a normalized vector representing respective signed amplitudes of their oscillations. The second mode shape and second natural frequency occur when the two masses are released after being displaced (by the correct distances) in opposite directions. In this case, the values in the mode shapes will have different signs, because the masses are moving in opposite directions.

Task #1
Write a single function script to perform the Third Order Adams-Bashforth-Moulton algorithm on a system of ordinary differential equations (IVP). The program must be able to handle any number of differential equations (an IVP with any number of state variables).
Your function m-file will receive through its header the following problem parameters (in this
exact order):
a. a cell array containing the function handles of the equations describing the time derivatives of the several state variables
b. the initial augmented state-array – a numerical row array containing the needed starting values of the independent variable and the several state variables (in the same order as the function handles in the cell array from the first input argument)
c. the step size to be used, and
d. the number of steps to take

Function Output
Your file should not return anything to the calling program. Instead, it should
i. print to screen (the command window) the final state of the system, including t, and
ii. plot a trace for each of the dependent variables as a function of time.
Plot each trace on a separate graph. To facilitate plotting, store the results of all time steps in the augmented state-array
a. each row represents the system at a different point in time
b. column 1 contains the time, and each other column contains the temporal history of a particular state variable
Store only the final result of each time step, not the predicted or non-converged corrected values. Title each plot ‘Var1’, ‘Var2’, etc. Try help plot and help title in the MATLAB command window to learn more about plotting. num2str will also be useful for creating the labels themselves.

Hints for your Function Script
Handling a System of Un-predetermined Size
Technical
Because you will not know in advance the order of the IVP your function will be asked to solve, some of the system parameters will need to be passed in as arrays, rather than as scalars. The step size and number of steps will still be scalars, but there will be an un-predetermined number of state variables, each of which will have an associated starting value and derivative expression. The starting value of each state variable is a scalar, so the collection of them can be stored in a normal array. On the other hand, your function script must be passed a separate function handle for the equation for the time derivative of each state variable. To store function handles, we must use a cell array. So then, you must be able to evaluate a function handle stored in a cell array. We did this in Lab #4, but I cover it again later in this document.
NOTE: Because your augmented state array will include time as the first value, the function stored in the nth index of the cell array is actually the derivative of the state variable in the (n+1)th index of the augmented state array.

Metaphysical
To figure out how to structure your code to deal with any size system, consult HW#9 and the notes from lab. In the HW, we performed RK3 on a 2nd order system. We can call the state variables x and v, and we can call their slopes (time derivatives) k and m, respectively. Feel free to write out and logically group the steps you took to solve the system for one timestep. Now, imagine we have 3 or 4 state variables each with their own slopes. Again, write out and logically group the steps you took to solve the system of equations. What changed about the process? Did the number of groups increase? Did the number of steps within particular groups increase? Now, we need to replace the RK4 algorithm steps with the ABM4 algorithm steps. Pay careful attention to the following differences:
a. the predictor and its modifier are performed once per time step, the corrector and its modifier are repeated until convergence each time step
b. the is no modifier on the predictor in the first time step
c. when performing the modifier steps, you use only non-modified values
d. when you check for convergence you compare only modified values

Function handles with unknown number of inputs
If you remember doing this for Modified Secant on a system of equations, then you may skip this section.
Up to now, you’ve probably only used function handles for equations that had a predetermined number of inputs. For this program, each derivative equation may be a function of the independent variable and any or all of the state variables. This is only a problem because you don’t know now (as you’re writing the function script) how many or which of the state variables will appear in each of the derivative equations. Fortunately, MATLAB allows us to put arrays into the equation function evaluation.
For example, if I define a function handle:
d1_by_dt = @(t,a) a(1) * a(3) / (2 + t);
when you ask MATLAB to evaluate the mathematical expression, MATLAB will expect to receive a scalar as the 1st argument, and an array of at least three terms as the 2nd argument. It will use the first and third terms in the array, along with the scalar, to evaluate the expression.
The augmented state array passed into the function already has t as the first value, so let’s just stick with that. Then our example equation becomes:
d1_by_dt = @(a) a(2) * a(4) / (2 + a(1));
You, of course, will have no prior knowledge of which state variables are used in each derivative expression. Therefore, for each derivative evaluation, you will need to submit an array containing the correct values of the dependent variable and all of the state variables, in the same order as the initial values in the passed-in augmented state array. When grading your code, it will be my responsibility to ensure that all of the derivative equations use the correct variables.

Function handles in cell arrays
The function handles of the derivative expressions of the several state variables will be passed to your function script in a cell array. The syntax for evaluating a function handle that is stored in a cell array is as follows:
value = array_name{n}(input1, input2, ... , inputQ)
• array_name is the name of the cell array in which the handles are stored
• n is the index into the cell array, in which is stored the expression for (d/dt)y n
• input1 thru inputQ are the inputs to the mathematical expression represented by the function handle – in our case that will simply be the augmented state array
Please note that those are braces around the index ‘n’ and parentheses around the arguments.

The ABM4 process
To calculate the predicted value of each state variable, you are going to perform very similar math:
n+1 ^0y k = n y k + h((55/24)g^k n - (59/24)g^k n-1 + (37/24)g^k n-2 - (3/8)g^k n-3)

where the k superscript on each g term matches the state variable number (the k subscript on each y term). Take advantage of how MATALB deals with vectors to simplify the process of calculating the predictor for all the state variable and each subsequent corrector for all the state variables. Create an array that holds the g m evaluations for all of the state variables. You will have 3 or 4 of these, one for each timestep. Then you can ...
temp_state = current_state + (h / 2) * slopes_array

Test your program
The spring-mass system from the Terminology section has been reproduced here.
Neglecting friction, the equations of motion of this system are:
ẍ1 = (-(k1+k2)x1 + k2x2)/m1
ẍ2 = (-(k2+k3)x2 + k2x1)/m2

where x1 and x2 are measured from the static equilibrium locations of the respective masses. To handle the 2nd order nature of the time derivatives, we must include additional state variables v1 and v2. If we arrange the state variables in the order x1, v1, x2, v2, the augmented state array becomes t, x1, v1, x2, v2, and the corresponding differential equations becomes:
dx1/dt = g1(a) = a(3)
dv1/dt = g2(a) = (-(k1+k2)a(2) + k2a(4))/m1
dx2/dt = g3(a) = a(5)
dv2/dt = g4(a) = (-(k2+k3)a(4) + k2a(2))/m2

Remember, a(1) is time, which does not appear in our equations of motion.
Create a main script (program) that performs the following tasks:
iii. define the system parameters (masses and spring stiffnesses)
iv. define the IVP functions above and store them in a cell array
v. set the initial condition of the augmented state array (time and all four state variables)
vi. call your RK4 function script
Set the system parameters to m1 = m2 = 1(kg) and k1 = k2 = k3 = 13.16(N/m). Set the initial
conditions to t = 0, x1(0) = 0.01(m), v1(0) = 0, x2(0) = -0.01(m), and v2(0) = 0. Call your RK4 function script with a time step of 0.001 seconds and 2000 time steps. Inspect the plots of x1 and x2. They should be mirrored sinusoidal functions with amplitudes of 0.01 and frequencies of 1Hz (2π rad/sec).

Analyze Your Custom Vibration System
You are going to customize the parameters of the spring mass system based upon your student ID number. Consider digits 4-5 as a two digit number X. Let digits 6-7 be a two digit number Y. If either X or Y are zero, pretend that they are 35.

Task #2: Apply Your Program
You will use your function to analyze a customized version of the vibration system depicted on the previous page. The system will be customized based on the last four digits of your student ID. Your goal will be to determine the second natural frequency and corresponding mode shape. For more details, refer to the Terminology section of this document.
Set the system parameters to m1 = 1(kg), m2 = (0.75-(X/200))(kg), k1 = k2 = 13.16(N/m), and k3 = (13.66+(Y/200))(N/m). Set the initial conditions to t = 0, x1(0) = -0.01(m), v1(0) = 0, and v2(0) = 0. Please note that the initial value of x1 is negative. Now, find the second natural frequency and second mode shape of your system by running your RK4 function multiple times, with different initial values for x2. I think the correct value should be > 0.01.
When you find the correct value of x2(0) each mass should demonstrate smooth sinusoidal motion. Inspect your table of state variables to determine the period of the oscillation. The period is the time between consecutive maximum displacements. The frequency in Hz is the inverse of the period. The mode shape is the vector [x1(0); x2(0)] scaled to have a Euclidean norm of unity (1).
Include in your Lab Report the following:
• the last four digits of your student ID number
• the values you used for X and Y
• the second natural frequency, in Hz, of your customized spring-mass system
• the normalized second mode shape of your customized spring-mass system
• plots of x1 and x2 as functions of time for two seconds
Provide the natural frequency and mode shape values to 5 significant figures.

Here is the RK4 chunk that is supposed to go on the top of the code for the function:
% In the code below:
%
% g is cell array containing the various dy/dt functions
%
% s_a is the augmented state array. It starts out as a row vector. It is
% passed in with the starting values of t and each state variable, in the
% same order as the function handles in the g cell array.  It will
% eventually contain all of the values of t and the state variables over
% time.
%
% h is the step size

N = length(g);
slopes = zeros(4,N)

for iter = 1:2
    state = s_a(iter,:);
    for c = 1:N
        slopes(1,c) = g{c}(state);
    end
    
    state = s_a(iter,:) + (h/2)*[1, slopes(1,:)];
    for c = 1:N
        slopes(2,c) = g{c}(state);
    end
    
    state = s_a(iter,:) + (h/2)*[1, slopes(2,:)];
    for c = 1:N
        slopes(3,c) = g{c}(state);
    end
    
    state = s_a(iter,:) + h*[1, slopes(3,:)];
    for c = 1:N
        slopes(4,c) = g{c}(state);
    end
    
    RK = (slopes(1,:) + 2*slopes(2,:) + 2*slopes(3,:) + slopes(4,:)) / 6;
    
    s_a(iter+1,:) = s_a(iter,:) + h*[1, RK]
        
end

When this chunk is done, we should have three rows in the augmented state array, corresponding to our first three time points. It shouldn't be changed in any way while making the function script.