Fit experimental data to analytical expression

Hi,

I'm struggling to fit an experimentally measured variable to an analytical expression. The experimentally measured variable is called the phase velocity (Eq 2; screenshot attached) and I am aiming to fit it within an analytical expression (Eq 3; also attached as a screenshot) to determine mu and eta, which are the shear modulus and shear viscosity of a sample. Now the paper I am using as reference describes doing the following: "The final phase velocity curve was then fitted by the linear least squares method to (2) and (3) to obtain the shear viscoelasticity of the tissues. The shear elasticity and shear viscosity were varied in the model y(mu, eta) with step size of 0.1 kPa and 0.1 Pa · s until the error between the model and the phase velocity data was minimized". I am struggling to implement this in MATLAB. My questions are the following:

1/ What is the best way to fit the experimentally measured phase velocity to Eq 3? Is it with symbolic variables or fminsearch or neither? How do I even write out equation 3 given that it is a two-sided equation and phase velocity is a function of variables on both sides? I have tried the following based on what I've read online (but I am sure it is incorrect):

k_L = omega ./ phase_vel;
numPar = 2;
mu = sym('mu', [1 numPar]);
thickness = 4e-3;
ks = omega .* sqrt(rho ./ (mu));
kl = omega ./ phase_vel;
h = thickness / 2;
beta = sqrt (kl.^2 - ks.^2);
first_term = 4 .* kl.^3 .* beta .* cosh ( kl .* h) .* sinh ( beta .* h);
second_term = ( ks.^2 - 2 .* ( kl.^2 )).^2 .* sinh ( kl .* h ).* cosh( beta .* h );
third_term = ks .^ 4 .* cosh ( kl .* h) .* cosh ( beta .* h);
fun1 = symfun(first_term - second_term - third_term, mu);
myfun = matlabFunction(fun1, 'Vars', {mu});
mu = [10e3, i .* omega .* 0.01];
[x, fval, exitflag, output] = fminsearch (mufun, mu);

Here, phase_vel would be the experimentally measured values.

2/ How would I determine the mu and eta from this using LLS to obtain the values? Is there a way to integrate it from above?

Thank you in advance for your time and help! Any code that would help me address 1/ and 2/ would be very appreciated. Have a great day! :)

[Cris LaPierre]

I would use the curve fitting app. Set the fit to Custom and enter your fit equation.

[Torsten]

You say you measured the phase velocity - so phase velocity seems to be the dependent variable in your model. What is the independent variable - thus which parameter was varied to get different phase velocities ? Is it omega ?

[Joseph]

I will look at this also but I'd to come up with an automated solution for multiple data files - wouldn't I need to run the app with each file? If so, wouldn't this be time-consuming?

Answer 1

[Star Strider]

Equation (3) is a bit mystifying to me.

imshow(imread('analytical_expr.png'))

imshow(imread('phase_vel.png'))

I assume that ‘phase_vel’ is the independent variable. (I can’t determine what the dependent variable is, or if it should be ‘phase_vel’ and something else should be the independent variable.)

The matlabFunction call needs to be:

myfun = matlabFunction(fun1, 'Vars', {[mu],phase_vel,omega,rho});

with the appropriate independent variable in the argument list as determined in the 'Vars' value (the square brackets around ‘mu’ will return it as a vector, probably called ‘in1’ that by convention must be a row vector), since in MATLAB the objective function has two arguments, the first being the parameter vector to be estimated, and the second the independent variable. Other additional parameters can be included in the argument list.

The fminsearch call is then:

[x, fval, exitflag, output] = fminsearch (@(mu)norm(dependent_variable - myfun(mu,phase_vel,omega,rho)), mu0);

where ‘mu0’ is:

mu0 = [10e3, i .* omega .* 0.01];

Clarification would be helpful.

[Joseph]

Thanks for your answer and apologies for the confusion. phase velocity is dependent on omega (angular frequency). Everything is angular-frequency-dependent (or omega-dependent). I'm not sure if the code I wrote is the best way to do this so if you know of a better way also, I'm open to that.

[Torsten]

If omega is the independent and phase velocity is the dependent variable, why is equation (2) not sufficient for fitting mu and eta ? Is equation (3) a kind of constraint on mu and eta that has to be satisfied ?

[Star Strider]

No apologies necessary! I appreciate the clarification.

In that instance, the matlabFunction call becomes:

myfun = matlabFunction(fun1, 'Vars', {[mu],omega,rho});

and the fminsearch call becomes:

[x, fval, exitflag, output] = fminsearch (@(mu)norm(phase_vel - myfun(mu,omega,rho)), mu0);

assuming that ‘phase_vel’ and ‘omega’ are equal-sized (preferably column) vectors already present in the workspace.

You will need to test this with the actual data.

[Joseph]

@Torsten In elastic materials, this is what is done. But in thin, linear, viscoelastic materials, Eq. 3 is used because it accounts for different material properties and geometries.

@Star Strider thank you! I will try this out and let you know if I have any questions. I really appreciate your quick replies and help. Thank you so much!

[Torsten]

So you want to minimize the sum of least squares simultaneously for equations (2) and (3) ?

Then you will have to use lsqnonlin since equation(3) cannot be solved explicitly for the dependent variable (phase velocity).

[Joseph]

Okay! Thank you.

[Star Strider]

@Joseph — My pleasure!

@Torsten — Noted.

[Joseph]

Would either of you be able to take a look at another question I asked? https://www.mathworks.com/matlabcentral/answers/2088211-taking-the-2d-fft-of-axial-displacement-data?s_tid=srchtitle

[Cris LaPierre]

You haven't accepted answers to any of the questions you have asked. If your questions have been answered, please accept the answer. It helps others know that question has a solution, and contributes to the reputation of the person who wrote the answer.

[Joseph]

Sorry about that Cris! I accepted the answer now.

[Joseph]

I accepted all the other answers also. Thanks for the reminder, Cris. Apologies again.

Any thoughts on my other question linked above?