EEU22E11 Curve fitting -Solved

The compressive properties of human muscle tissue are significant for impact biomechanics, surgical simulation and rehabilitation engineering. The most important loading direction is transverse to the muscle fibre direction, and tests performed at TCD on pig muscle samples using a universal testing machine showed nonlinear anisotropic behavior under quasi-static loading conditions 1. In this assignment you will fit regression

curves to data from a compression test on a sample performed in the transverse (or cross fibre XF) direction, i.e. corresponding to the blue curve in the figure below. You will see that the relationship between stress and strain is not linear, and therefore defining a constant value for slope of the stress strain curve (ie assuming that there is a Young’s Modulus) will not give a very good fit to the data.

The data is stored in a file muscle_data_2017.xls. Use the following Matlab code to read in the data.

raw_data = xlsread(’muscle_data_2017.xlsx’); strain = raw_data(:,3); stress = raw_data(:,4);

You will now explore two different models for the data. We will use the model formulation yˆk = f(φk) where yˆk is our estimated stress given the measured strain, φk, for the kth measurement. The error between the observed stress yk and the modelled (estimated) stress yˆk is therefore ek = yk−yˆk.

1 M. Van Loocke, C.G. Lyons, and C.K. Simms. A validated model of passive muscle in compression. Journal of Biomechanics, pages 2999–3009, October 2005

1.   Plot in figure 1, a plot of the observed stress (on the yaxis) yk Vs strain (on the x-axis) (φk). Use a blue line of thickness 3.

2.   You are required to fit the data to a polynomial model. The problem is to estimate the best or appropriate model order m. Test 5 model orders m = [1 : 5]. For each polynomial model with order m use the polyfit function in Matlab to estimate the coefficients of the polynomial fit to the data. Calculate the sum squared error for each polynomial fit and assign the resulting data to the vector sse_per_m. That vector should be 5 elements long because you will have tested 5 model orders. In figure 2, plot sse_per_m versus model order and use that to select the optimal model order mˆ for this problem. Use a red line of thickness 3. Label the figure and axes appropriately. Assign your selected model order to the variable my_m.

3.   Using your selected model order mˆ , calculate the estimated signal yˆk and superimpose that line in red on the plot in figure 1. Assign your estimated signal yˆk to the variable est_stress.

You will have to use your notes or the

 
book by Chapra to find the definition
4. Calculate the error between your estimate and the observed stress, plot a histogram in figure 3 showing the distribution of errors. Your plot should look something like the plot
for the coefficient of determination.
shown on the side panel (no need for the title though). Note
Std error = 0.010037r = 0.99943
 that the values shown in that plot may be different from your estimates. Calculate the coefficient of determination for this polynomial fit and assign that to the variable ccoef_p.

5.    Now fit the data to a different nonlinear model as follows.

(You will have to use fminsearch see the notes below).

                                           yˆk = 1 −a2− (a0 exp(a1φk))                        (1)                      Error in Polynomial fit

 where yˆ is the estimated stress given a strain measurement φ. You must start from an initial estimate [0.3,0.3,0.3]. Using your estimated values for a0,a1,a2 calculate the estimated signal yˆk and superimpose that line in black/dashed on the plot in figure 1. Your final figure 1 should look something like the plot in the side panel shown here.

6.    Calculate the error between your estimate and the observed stress, and assign it to the variable err_nl . Plot a histogram in figure 4 showing the distribution of errors. (Just use the hist command in Matlab e.g. hist(err_nl)).

Hence calculate the coefficient of determination for this nonlinear fit, and assign that to the variable ccoef_nl.

Notes
Use » help fminsearch and » help polyfit from the matlab command line (») to find out more about these matlab functions and how to use them. For fminsearch you will need to pass arguments to your optimisation function which are not being estimated i.e. the data itself stress, strain. You may do this using the syntax coefs = fminsearch(@(x) func(x, stress, strain),[0.3;0.3;0.3]);. In this invocation you are declaring to fminsearch that the error measurement is in the function func(), and the @(x) is declaring that only the variable x in this function is being optimised. You will have to define the function func() yourself. It will be based on the model you are testing. In this case x should be a vector with your three model parameters x(1) = a0,x(2) = a1,x(3) = a2. Hence func() then

should calculate the sum squared error between the estimated stress using the model given (using your parameters x) and the actual measured stress, given the measured strain.

You might find it useful to know that Matlab string handling uses vector syntax like this.

str = [’This variable =’, num2str(number), ’while that = ’, num2str(value)];

The num2str function turns a number into a string for use in string handling like this. The string itself is in str and it just concatenates all the strings in the vector. You can then use this string for auto-composition of a title. For instance e.g. title(str, ’fontsize’ , 20) will create a title that has the text as follows.

This variable = 78.2 while that = 7.1