Fitting constants of a matrix to experimental data

Hello, I am dealing with matrices where the initial matrix is [lam1, 0, 0; 0, lam2, 0; 0, 0, lam3]

The function is as follows:

I have experimental data and i would like to fit parameters G, K, k1, k2, gamma, to the data but solve lam2 and lam3 such that sigma_cauchy(2,2) and sigma_cauchy(3,3) are set to zero and just vary lam1. Any ideas how to set this up would be greatly appreciated thank you?

D_1 = 2/K

function sigma_cauchy_fun_nom = sigma_cauchy_fun(G, K, k1, k2, gamma, lam2, lam3, lam1)

    F = [lam1, 0, 0; 0 lam2, 0; 0, 0, lam3];
    J = det(F); %determinant of deformation gradient
    C = transpose(F)*F; %right cauchy green tensor
    B = F*transpose(F); %left cauchy green tensor
    I_1 = trace(B); %Inveriant
    sigma_iso = (G/(J^(5/3))) * (B - ((I_1/3)*I)) + (K*(J-1)*I)
    sigma_vol = (2/D_1) * (J - 1) * I;
    sigma_iso_vol = sigma_iso + sigma_vol;
    m_4 = [cosd(gamma); sind(gamma); 0];
    m_6 = [cosd(-gamma); sind(-gamma); 0];
    I_4 = transpose(m_4) * C * m_4;
    I_6 = transpose(m_6) * C * m_6;
    a_4 = F*m_4;
    a_6 = F*m_6;
    sigma_aniso = (2*k1*(I_4 - 1)*exp(k2*((I_4 - 1)^2))*(a_4*transpose(a_4))) + (2*k1*(I_6 - 1)*exp(k2*((I_6 - 1)^2))*(a_6*transpose(a_6)));
    sigma_cauchy = sigma_iso_vol + sigma_aniso;
    P11 = J * sigma_cauchy * inv(transpose(F));
    P_nom_stress(i) = P11(1,1);

end

[Torsten]

Unclear.

[Oisín Conway]

@Torsten For example at the moment I am guessing parameters G, K, k1, k2, gamma but how would you use fittype() and fit() using a custom defined function like below to solve for the parameters G, K, k1, k2, gamma? thanks

for i = 1:np
    
    P11_tensor = sigma_cauchy_fun(G, K, k1, k2, gamma, lam1(i))
    
    P_nom_stress(i) = P11_tensor(1,1)
    
    
end 

figure(1)
clf
plot(strain_axial_exp_nom, stress_axial_exp_nom, 'color', 'k')
xlabel('Strain (mm/mmm)')
ylabel('Stress (MPa)')
xlim([0 max([max(strain_axial_exp_nom), max(strain_circum_exp_nom)])])
hold on
%plot(strain_circum_exp_nom, stress_circum_exp_nom, 'color', 'b')
plot(lam1-1, P_nom_stress)
%plot(lam1-1, P_nom_stress_circum, 'r')
legend("Axial (Expermential)", "Circumferential (Experimental)", "Axial (Theoretical)", "Circumferential (Theoretical)", 'location', 'northwest')
hold off

function P11 = sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)
    
    I = eye(3,3);
    D_1 = 2/K;
    F = [lam1, 0, 0; 0 1/sqrt(lam1), 0; 0, 0, 1/sqrt(lam1)];
    J = det(F); %determinant of deformation gradient
    C = transpose(F)*F; %right cauchy green tensor
    B = F*transpose(F); %left cauchy green tensor
    I_1 = trace(B); %Inveriant
    sigma_iso = (G/(J^(5/3))) * (B - ((I_1/3)*I)) + (K*(J-1)*I);
    sigma_vol = (2/D_1) * (J - 1) * I;
    sigma_iso_vol = sigma_iso + sigma_vol;
    m_4 = [cosd(gamma); sind(gamma); 0];
    m_6 = [cosd(-gamma); sind(-gamma); 0];
    I_4 = transpose(m_4) * C * m_4;
    I_6 = transpose(m_6) * C * m_6;
    a_4 = F*m_4;
    a_6 = F*m_6;
    sigma_aniso = (2*k1*(I_4 - 1)*exp(k2*((I_4 - 1)^2))*(a_4*transpose(a_4))) + (2*k1*(I_6 - 1)*exp(k2*((I_6 - 1)^2))*(a_6*transpose(a_6)));
    sigma_cauchy = sigma_iso_vol + sigma_aniso;
    P11 = J * sigma_cauchy * inv(transpose(F));
%     P_nom_stress(i) = P11(1,1);

end

[Torsten]

See the example

Fit a Curve Defined by a File

under

https://de.mathworks.com/help/curvefit/fit.html

Or a simple one here:

x = (0:0.1:2).';
y = x.^2 + 0.01*(-1+2*rand(size(x)));
fun = @(a,b,c,x) fitting(a,b,c,x);
ft = fittype(fun);
f = fit( x, y, ft, 'StartPoint',[1 1 1])

f = 
     General model:
     f(x) = fitting(a,b,c,x)
     Coefficients (with 95% confidence bounds):
       a =   0.0002479  (-0.007652, 0.008148)
       b =   -0.001645  (-0.01995, 0.01666)
       c =           1  (0.9915, 1.009)

plot(f,x,y)

function y = fitting(a,b,c,x)
  y = a + b*x + c*x.^2;
end

[Oisín Conway]

Thank you so much for that example @Torsten would this work for coefficients of a matrix as well? for example:

where G, K, k1, k2 and gamma coefficients and lam1 is independent

function P_nom_stress = sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)
    
    I = eye(3,3);
    D_1 = 2/K;
    F = [lam1, 0, 0; 0 1/sqrt(lam1), 0; 0, 0, 1/sqrt(lam1)];
    J = det(F); %determinant of deformation gradient
    C = transpose(F)*F; %right cauchy green tensor
    B = F*transpose(F); %left cauchy green tensor
    I_1 = trace(B); %Inveriant
    sigma_iso = (G/(J^(5/3))) * (B - ((I_1/3)*I)) + (K*(J-1)*I);
    sigma_vol = (2/D_1) * (J - 1) * I;
    sigma_iso_vol = sigma_iso + sigma_vol;
    m_4 = [cosd(gamma); sind(gamma); 0];
    m_6 = [cosd(-gamma); sind(-gamma); 0];
    I_4 = transpose(m_4) * C * m_4;
    I_6 = transpose(m_6) * C * m_6;
    a_4 = F*m_4;
    a_6 = F*m_6;
    sigma_aniso = (2*k1*(I_4 - 1)*exp(k2*((I_4 - 1)^2))*(a_4*transpose(a_4))) + (2*k1*(I_6 - 1)*exp(k2*((I_6 - 1)^2))*(a_6*transpose(a_6)));
    sigma_cauchy = sigma_iso_vol + sigma_aniso;
    P11 = J * sigma_cauchy * inv(transpose(F));
    P_nom_stress = P11(1,1);

end

[Torsten]

would this work for coefficients of a matrix as well?

Why not ? But you had to loop because lam1 will most probably be a vector as input and you had to return a vector for P_nom_stress as output.

[Oisín Conway]

I have it now as follows (not in a loop), strain_axial_exp_nom and stress_axial_exp_nom are column vectors. but it throws an error:

Error using fittype>iAssertIsMember (line 1084)

Coefficient gamma does not appear in the equation expression.

Error in fittype>iTestCustomModelParameters (line 731)

iAssertIsMember( obj.coeff, variables, 'curvefit:fittype:noCoeff' );

Error in fittype>iCreateFittype (line 372)

iTestCustomModelParameters( obj );

Error in fittype (line 330)

obj = iCreateFittype( obj, varargin{:} );

Error in bio_function_file_axial (line 80)

ft = fittype(fun, 'indep','lam1')

fun = @(G, K, k1, k2, gamma, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)
ft = fittype(fun, 'indep','lam1')
%ft = fittype(@(G, K, k1, k2, gamma, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1),'indep','lam1 ') 
f = fit(strain_axial_exp_nom, stress_axial_exp_nom, ft, 'StartPoint', [0 0 0 0 0])

[Torsten]

Not possible this way. Both parameter lists must coincide:

fun = @(G, K, k1, k2, gamma, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)

[Oisín Conway]

What do mean by this sorry? the parameters are exactly the same.

[Torsten]

A moment ago, gamma was missing - and that's what the error message also confirms:

fun = @(G, K, k1, k2, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)

Error using fittype>iAssertIsMember (line 1084)

Coefficient gamma does not appear in the equation expression.

[Oisín Conway]

Apologies, i realised that and ran it again. At the moment i have the following below. I am unsure how to correct this.

The errors i am getting now are:

Error using fittype/testCustomModelEvaluation (line 12)

Expression sigma_cauchy_fun(G,K,k1,k2,gamma,lam1) is not a valid MATLAB expression, has non-scalar coefficients, or cannot be evaluated:

Error while trying to evaluate FITTYPE function :

Dimensions of arrays being concatenated are not consistent.

Error in fittype>iCreateFittype (line 373)

testCustomModelEvaluation( obj );

Error in fittype (line 330)

obj = iCreateFittype( obj, varargin{:} );

Error in bio_function_file_axial (line 72)

ft = fittype(fun, 'independent', 'lam1', 'coefficients', {'G','K', 'k1', 'k2', 'gamma'})

Caused by:

Error using fittype/evaluate (line 88)

Error while trying to evaluate FITTYPE function :

Dimensions of arrays being concatenated are not consistent.

fun = @(G, K, k1, k2, gamma, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1)
ft = fittype(fun, 'independent', 'lam1', 'coefficients', {'G','K', 'k1', 'k2', 'gamma'})
%ft = fittype(@(G, K, k1, k2, gamma, lam1) sigma_cauchy_fun(G, K, k1, k2, gamma, lam1),'indep','lam1 ') 
f = fit(strain_axial_exp_nom, stress_axial_exp_nom, ft, 'StartPoint', [0 0 0 0 0])

[Torsten]

And did you loop over lam1 which is a vector input from fit to "sigma_cauchy_fun" ? And did you return a column vector for "P_nom_stress" ?

[Oisín Conway]

No i took it out of the for loop because it was not working.

[Torsten]

The input of fit to "sigma_cauchy_fun" will be a column vector of values for lam1 (which you must supply in the call to fit) and expects a vector of the same size as output (P_nom_stress) .

So your call to "fit" as

f = fit(strain_axial_exp_nom, stress_axial_exp_nom, ft, 'StartPoint', [0 0 0 0 0])

makes no sense - the first argument must be lam1 and the second argument the data vector you want to fit with the function generating the vector "P_nom_stress" (usually a measurement vector).

[Oisín Conway]

I see, so should it be the following then

lam1 = linspace(1, 1 + 0.6, 1212)
lam1 = transpose(lam1)
f = fit(lam1 - 1, stress_axial_exp_nom, ft, 'StartPoint', [0 0 0 0 0])