Inverse Problem, I have unknowns to solve it

I can make as many as equation I want with changing Kf, which it will give me A, B, C, D, F equal to the length of Kf, menas you can make it overdetermind

I want to find the values of A1,A2,A3,A4,A5, L and N

clc;
clear;
close all;
K1 = 60;
K2 = 40;
U1 = 30;
U2 = 20;
pt1 = 0;
pt2 = 0.4;
Kf = linspace(1, 3, 3);
Kf1 = 0;
Xl1 = 0.5;
Xl2 = 1-Xl1;
Uf1 = 0;
Uf2 = 0;
phi = Xl1 * pt1 + Xl2 * pt2;
[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);
[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);
Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;
Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;
[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);
Sij_A_Frame = inv(Cij_A_Frame);
Sig_D  = zeros(6, 6, numel(Kf));
for i = 1:numel(Kf)
    Kf2 = Kf(i);
    [~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);
    [~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);
    Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;
    Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;
    [Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);
    Sij_star = inv(Cij_A_Ud);
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;
end
A = squeeze(Sig_D(1,1,:)).';
B = squeeze(Sig_D(1,1,:)).';
C = squeeze(Sig_D(3,3,:)).'; 
D = squeeze(Sig_D(1,3,:)).';
F = squeeze(Sig_D(4,4,:)).';
kappa_f = 1./Kf;
kappa_a = trace(Sij_A_Frame);

%FinD A1,A2,A3,A4,A5, L and N
A = A1/((kappa_f - L) * phi + (kappa_a - N));

Unrecognized function or variable 'A1'.

B = A2/((kappa_f - L) * phi + (kappa_a - N));
C = A3/((kappa_f - L) * phi + (kappa_a - N));
D = A4/((kappa_f - L) * phi + (kappa_a - N));
F = A5/((kappa_f - L) * phi + (kappa_a - N));

function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)
Umin=min(U);
Umax=max(U);
Kmin=min(K);
Kmax=max(K);
if P==0
    Klb = 1./((1-P).*sum(Vfr./K));
elseif P>0
    Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));
end
Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;
Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));
if P ==0
    Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));
elseif P>0
    Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));
end
Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;
Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;
end 
function [c]=Backus_average(VF,G,Y)
x = sum(VF.*G.*(Y+G)./(Y+2*G));
y = sum(VF.*G.*Y./(Y+2*G));
z = sum(VF.*Y./(Y+2*G));
u = sum(VF./(Y+2*G));
v = sum(VF./G);
w = sum(VF.*G);
A = 4*x + z*z/u;
B = 2*y + z*z/u;
E = 1/u;
F = z/u;
D = 1/v;
M = w;
c = zeros(6,6);
c(1,1) = A;
c(1,2) = B;
c(1,3) = F;
c(2,1) = B;
c(2,2) = A;
c(2,3) = F;
c(3,1) = F;
c(3,2) = F;
c(3,3) = E;
c(4,4) = D;
c(5,5) = D;
c(6,6) = M;
end

Answer 1

[Torsten]

L and N cannot be estimated separately.

As long as phi is a single scalar value , L*phi + N must be treated as only one parameter that replaces L*phi + N. Otherwise, you do overfitting.

I didn't do that in the code below, but it should be easy for you to modify it accordingly.

clc;
clear;
close all;
K1 = 60;
K2 = 40;
U1 = 30;
U2 = 20;
pt1 = 0;
pt2 = 0.4;
Kf = linspace(1, 3, 3);
Kf1 = 0;
Xl1 = 0.5;
Xl2 = 1-Xl1;
Uf1 = 0;
Uf2 = 0;
phi = Xl1 * pt1 + Xl2 * pt2;
[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);
[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);
Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;
Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;
[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);
Sij_A_Frame = inv(Cij_A_Frame);
Sig_D  = zeros(6, 6, numel(Kf));
for i = 1:numel(Kf)
    Kf2 = Kf(i);
    [~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);
    [~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);
    Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;
    Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;
    [Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);
    Sij_star = inv(Cij_A_Ud);
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;
end
A = squeeze(Sig_D(1,1,:)).';
B = squeeze(Sig_D(1,1,:)).';
C = squeeze(Sig_D(3,3,:)).'; 
D = squeeze(Sig_D(1,3,:)).';
F = squeeze(Sig_D(4,4,:)).';
kappa_f = 1./Kf;
kappa_a = trace(Sij_A_Frame);
x = [A,B,C,D,F];
y = @(A1,A2,A3,A4,A5,L,N)[A1./((kappa_f - L) * phi + (kappa_a - N)),...
                          A2./((kappa_f - L) * phi + (kappa_a - N)),...
                          A3./((kappa_f - L) * phi + (kappa_a - N)),...
                          A4./((kappa_f - L) * phi + (kappa_a - N)),...
                          A5./((kappa_f - L) * phi + (kappa_a - N))];
fy = @(p)y(p(1),p(2),p(3),p(4),p(5),p(6),p(7));
F = @(p) fy(p) - x;
sol = lsqnonlin(F,[7 6 5 4 3 2 1])

Local minimum found.

Optimization completed because the size of the gradient is less than
the value of the optimality tolerance.

sol = 1×7

0.0360 0.0360 0.0981 -0.0222 0.3554 33.9358 8.5194

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%FinD A1,A2,A3,A4,A5, L and N
%A = A1/((kappa_f - L) * phi + (kappa_a - N));
%B = A2/((kappa_f - L) * phi + (kappa_a - N));
%C = A3/((kappa_f - L) * phi + (kappa_a - N));
%D = A4/((kappa_f - L) * phi + (kappa_a - N));
%F = A5/((kappa_f - L) * phi + (kappa_a - N));

function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)
Umin=min(U);
Umax=max(U);
Kmin=min(K);
Kmax=max(K);
if P==0
    Klb = 1./((1-P).*sum(Vfr./K));
elseif P>0
    Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));
end
Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;
Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));
if P ==0
    Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));
elseif P>0
    Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));
end
Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;
Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;
end 
function [c]=Backus_average(VF,G,Y)
x = sum(VF.*G.*(Y+G)./(Y+2*G));
y = sum(VF.*G.*Y./(Y+2*G));
z = sum(VF.*Y./(Y+2*G));
u = sum(VF./(Y+2*G));
v = sum(VF./G);
w = sum(VF.*G);
A = 4*x + z*z/u;
B = 2*y + z*z/u;
E = 1/u;
F = z/u;
D = 1/v;
M = w;
c = zeros(6,6);
c(1,1) = A;
c(1,2) = B;
c(1,3) = F;
c(2,1) = B;
c(2,2) = A;
c(2,3) = F;
c(3,1) = F;
c(3,2) = F;
c(3,3) = E;
c(4,4) = D;
c(5,5) = D;
c(6,6) = M;
end

[Abdolrazzagh]

Thank you very much, could you please edit that based on your suggestion so I can use this for other problms that I have,

Thank you again

[Torsten]

clc;
clear;
close all;
K1 = 60;
K2 = 40;
U1 = 30;
U2 = 20;
pt1 = 0;
pt2 = 0.4;
Kf = linspace(1, 3, 3);
Kf1 = 0;
Xl1 = 0.5;
Xl2 = 1-Xl1;
Uf1 = 0;
Uf2 = 0;
phi = Xl1 * pt1 + Xl2 * pt2;
[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);
[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);
Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;
Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;
[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);
Sij_A_Frame = inv(Cij_A_Frame);
Sig_D  = zeros(6, 6, numel(Kf));
for i = 1:numel(Kf)
    Kf2 = Kf(i);
    [~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);
    [~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);
    Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;
    Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;
    [Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);
    Sij_star = inv(Cij_A_Ud);
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;
end
A = squeeze(Sig_D(1,1,:)).';
B = squeeze(Sig_D(1,1,:)).';
C = squeeze(Sig_D(3,3,:)).'; 
D = squeeze(Sig_D(1,3,:)).';
F = squeeze(Sig_D(4,4,:)).';
kappa_f = 1./Kf;
kappa_a = trace(Sij_A_Frame);
x = [A,B,C,D,F];
y = @(A1,A2,A3,A4,A5,LN)[A1./(kappa_f  * phi + kappa_a - LN),...
                         A2./(kappa_f  * phi + kappa_a - LN),...
                         A3./(kappa_f  * phi + kappa_a - LN),...
                         A4./(kappa_f  * phi + kappa_a - LN),...
                         A5./(kappa_f  * phi + kappa_a - LN)];
fy = @(p)y(p(1),p(2),p(3),p(4),p(5),p(6));
F = @(p) fy(p) - x;
sol = lsqnonlin(F,[7 6 5 4 3 2])

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 1×6

0.0329 0.0327 0.0873 -0.0189 0.3137 13.5292

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

%FinD A1,A2,A3,A4,A5, L and N
%A = A1/((kappa_f - L) * phi + (kappa_a - N));
%B = A2/((kappa_f - L) * phi + (kappa_a - N));
%C = A3/((kappa_f - L) * phi + (kappa_a - N));
%D = A4/((kappa_f - L) * phi + (kappa_a - N));
%F = A5/((kappa_f - L) * phi + (kappa_a - N));

function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)
Umin=min(U);
Umax=max(U);
Kmin=min(K);
Kmax=max(K);
if P==0
    Klb = 1./((1-P).*sum(Vfr./K));
elseif P>0
    Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));
end
Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;
Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));
if P ==0
    Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));
elseif P>0
    Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));
end
Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;
Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;
end 
function [c]=Backus_average(VF,G,Y)
x = sum(VF.*G.*(Y+G)./(Y+2*G));
y = sum(VF.*G.*Y./(Y+2*G));
z = sum(VF.*Y./(Y+2*G));
u = sum(VF./(Y+2*G));
v = sum(VF./G);
w = sum(VF.*G);
A = 4*x + z*z/u;
B = 2*y + z*z/u;
E = 1/u;
F = z/u;
D = 1/v;
M = w;
c = zeros(6,6);
c(1,1) = A;
c(1,2) = B;
c(1,3) = F;
c(2,1) = B;
c(2,2) = A;
c(2,3) = F;
c(3,1) = F;
c(3,2) = F;
c(3,3) = E;
c(4,4) = D;
c(5,5) = D;
c(6,6) = M;
end

[Abdolrazzagh]

I am afarid to let you know this L*phi + N making one unkown is wrong as we are changing the equation in this case.

However, the first one is correct, Why this solution is highlt depended n the intial guess, is there anyway to make this to not dependent on the initial guess, Also why lsqnonlin, and not fmincon?

[Torsten]

Every combination of L and N that makes LN = L*phi + N solves the fitting task. You must consider L*phi + N as one constant as long as phi is a constant that is equally valid for A, B, C, D and F.

"lsqnonlin" is especially suited for least-squares problems - that's why I chose it. If you want to set constraints on the parameters, you will have to use "fmincon" instead.

Answer 2

[Abdolrazzagh]

What if we only solve

A1./((kappa_f - L) * phi + (kappa_a - N)​​

treating L and N as the unknowns with three equations? Once we have L and N, we can calculate the other variables.

Obviously, this approach should also work if we change the equation to

A2./((kappa_f - L) * phi + (kappa_a - N)

The calculated values of L and N should be similar in both cases or if I use the other cases.

but the problem they do not as I tried,

I have seem you are mathematician or wotks with mathematic. could you plese give your opinion and why it happens? ([e-mail]) if we can have an online meeting to clearify this problem and solve (would be appriciated)

Answer 3

[Walter Roberson]

A = squeeze(Sig_D(1,1,:)).';
B = squeeze(Sig_D(1,1,:)).';
C = squeeze(Sig_D(3,3,:)).'; 
D = squeeze(Sig_D(1,3,:)).';
F = squeeze(Sig_D(4,4,:)).';

each of those is a row vector with the same number of elements as Kf has

%FinD A1,A2,A3,A4,A5, L and N
A = A1/((kappa_f - L) * phi + (kappa_a - N));
B = A2/((kappa_f - L) * phi + (kappa_a - N));
C = A3/((kappa_f - L) * phi + (kappa_a - N));
D = A4/((kappa_f - L) * phi + (kappa_a - N));
F = A5/((kappa_f - L) * phi + (kappa_a - N));

In each case, A, B, C, D, F are known, but A1, A2, A3, A4, A5, L and N are unknown. These should instead be equations expressing that the left hand side is equal to the right hand side.

Notice though that the right hand sides all involve division by the same factor. So if you were to do something like

TEMP = ((kappa_f - L) * phi + (kappa_a - N));

then you would have

A == A1/TEMP
B == A2/TEMP
C == A3/TEMP
D == A4/TEMP
F == A5/TEMP

At this point we need to pause and ask whether the / operator is really the one intended, or if instead the ./ operator is intended, which in turn revolves around the expected size of TEMP and the expected size of A1 and so on.

We see that TEMP is constructed in part from kappa_f and kappa_f is constructed from 1./Kf and Kf is a row vector -- so kappa_f is a row vector and so TEMP must be a row vector. The A B etc. are all row vectors, so we have row_vector == UNKNOWN / row_vector . Experimentally, scalar / row_vector fails, and row_vector / row_vector yields a scalar rather than a row vector, so we deduce that the stated / operations must have been intended to be ./ operations. Unfortunately, we still cannot deduce whether the operation is row_vector == scalar ./ row_vector or if it is row_vector == row_vector ./ row_vector. However, if we rearrange to row_vector .* TEMP where TEMP is a row vector, then we can see that the result must be a row vector. This leads to

A1 = A .* TEMP;
A2 = B .* TEMP;
A3 = C .* TEMP;
A4 = D .* TEMP;
A5 = F .* TEMP;

Now we have a problem: TEMP is defined in terms of the underfined variables L and N, There are 5(*3) equations in 5(*3) unknowns excluding L and N... so there just are not equations to be able to deduce L or N. There is also not enough information to be able to deduce the expected size of L and N.

So... the system is not solvable.

[Abdolrazzagh]

Thanks.

[Torsten]

@Walter Roberson

As far as I understood the task, A, B, C, D and F can become arbitrary large. Thus we have a fitting problem where the number of unknows is 6 (or 7) and the number of equations is 5*N.

[Abdolrazzagh]

I made the problem in different way

may change the way to solve and make it better

what do you think?

this is the way that should be solved though, I made that eqautions from it as it was easier but it has a lot of ucertaintly

clc;

clear;

close all;

K1 = 60;

K2 = 40;

U1 = 30;

U2 = 20;

pt1 = 0;

pt2 = 0.4;

Kf = linspace(1, 7, 7);

Kf1 = 0;

Xl1 = 0.5;

Xl2 = 1-Xl1;

Uf1 = 0;

Uf2 = 0;

phi = Xl1 * pt1 + Xl2 * pt2;

[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);

[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);

Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;

Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;

[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);

Sij_A_Frame = inv(Cij_A_Frame);

Sig_D = zeros(6, 6, numel(Kf));

for i = 1:numel(Kf)

Kf2 = Kf(i);

[~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);

[~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);

Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;

Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;

[Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);

Sij_star = inv(Cij_A_Ud);

Sig_D(:,:,i) = Sij_A_Frame - Sij_star;

end

kappa_f = 1./Kf;

kappa_a = trace(Sij_A_Frame);

% solve this function for unknowns A, B, F, C, D, N, L

Sig = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L);

function Sig = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L)

S_A = zeros(6,6);

S_A(1,1) = A;

S_A(2,2) = A;

S_A(1,2) = B;

S_A(2,1) = B;

S_A(1,3) = F;

S_A(2,3) = F;

S_A(3,1) = F;

S_A(3,2) = F;

S_A(3,3) = C;

S_A(4,4) = D;

S_A(5,5) = D;

S_A(6,6) = (1/2) * (A-B);

Sig = S_A / ((kappa_f - L) * phi + (kappa_a - N));

end

function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)

Umin=min(U);

Umax=max(U);

Kmin=min(K);

Kmax=max(K);

if P==0

Klb = 1./((1-P).*sum(Vfr./K));

elseif P>0

Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));

end

Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;

Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));

if P ==0

Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));

elseif P>0

Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));

end

Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;

Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;

end

function [c]=Backus_average(VF,G,Y)

x = sum(VF.*G.*(Y+G)./(Y+2*G));

y = sum(VF.*G.*Y./(Y+2*G));

z = sum(VF.*Y./(Y+2*G));

u = sum(VF./(Y+2*G));

v = sum(VF./G);

w = sum(VF.*G);

A = 4*x + z*z/u;

B = 2*y + z*z/u;

E = 1/u;

F = z/u;

D = 1/v;

M = w;

c = zeros(6,6);

c(1,1) = A;

c(1,2) = B;

c(1,3) = F;

c(2,1) = B;

c(2,2) = A;

c(2,3) = F;

c(3,1) = F;

c(3,2) = F;

c(3,3) = E;

c(4,4) = D;

c(5,5) = D;

c(6,6) = M;

end

[Walter Roberson]

@Torsten

No,

A = squeeze(Sig_D(1,1,:)).';

and

for i = 1:numel(Kf)
%...
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;

so the third dimension of Sig_D is numel(Kf) so A and the rest are 1 x numel(Kf) rather than being something potentially large.

Meanwhile, the factor I noted as TEMP is 1x3 because it is formed from 1./Kf which is 1x3.

So we have 1x3 == A1/1x3 which only works if the / operator is really the ./ operator, in which case we have 1x3 == A1./1x3 and re-arranging, A1 == 1x3 .* 1x3 . There are 5 such equations, each with different (known) 1 x 3 on the left of the .*, using the same (not completely known) 1 x 3 on the right hand side. So we have 5*3 known values and 5*3 + numel(L) + numel(N) unknown values, which is not something you can fit, even if you restrict the unknown L and N to each be scalars.

[Abdolrazzagh]

I think I explained this incorrectly before, so let me restate it more clearly.

We compute Sig_D, which is a (6,6) matrix. This matrix is a function of Kf, meaning that changing Kf produces a different Sig_D. We can generate as many such matrices as we want by varying Kf; these matrices represent the measured data.

Next, we define the model prediction

Sig_D_P = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L);

where the parameters (A, B, C, D, F, N,L) are unknown.

to solve this we consider

The matrix Sig_D_P = Sig_D

Here, (S_A) is a (6, 6) matrix, while all parameters in the denominator are scalars.

Our goal is to determine:

• the five parameters defining the matrix (S_A): (A, B, C, D, F), and
• the two scalar parameters in the denominator: (L) and (N).

By varying Kf (and increasing the number of measurements), we obtain multiple measured matrices Sig_D. This allows us to minimize the error between Sig_D and Sig_D_P, enabling a more accurate estimation of the unknown parameters.

all the unknowns are constant although we vary the Kf and make a lot of Sig_D

clc;
clear;
close all;
K1 = 60;
K2 = 40;
U1 = 30;
U2 = 20;
pt1 = 0;
pt2 = 0.4;
Kf = linspace(1, 7, 7);
Kf1 = 0;
Xl1 = 0.5;
Xl2 = 1-Xl1;
Uf1 = 0;
Uf2 = 0;
phi = Xl1 * pt1 + Xl2 * pt2;
[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);
[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);
Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;
Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;
[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);
Sij_A_Frame = inv(Cij_A_Frame);
Sig_D  = zeros(6, 6, numel(Kf));
for i = 1:numel(Kf)
    Kf2 = Kf(i);
    
    [~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);
    [~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);
    
    Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;
    Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;
    
    [Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);
    Sij_star = inv(Cij_A_Ud);
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;
end
kappa_f = 1./Kf;
kappa_a = trace(Sij_A_Frame);
% solve this function for unknowns A, B, F, C, D, N, L
%Sig_D_P = Sig_D; 
%Sig_D_P = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L);
function Sig_D_P = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L)
S_A = zeros(6,6);
S_A(1,1) = A;
S_A(2,2) = A;
S_A(1,2) = B;
S_A(2,1) = B;
S_A(1,3) = F;
S_A(2,3) = F;
S_A(3,1) = F;
S_A(3,2) = F;
S_A(3,3) = C;
S_A(4,4) = D;
S_A(5,5) = D;
S_A(6,6) = (1/2) * (A-B);
Sig_D_P = S_A / ((kappa_f - L) * phi + (kappa_a - N));
end
function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)
Umin=min(U);
Umax=max(U);
Kmin=min(K);
Kmax=max(K);
if P==0
    Klb = 1./((1-P).*sum(Vfr./K));
elseif P>0
    Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));
end
Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;
Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));
if P ==0
    Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));
elseif P>0
    Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));
end
Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;
Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;
end 
function [c]=Backus_average(VF,G,Y)
x = sum(VF.*G.*(Y+G)./(Y+2*G));
y = sum(VF.*G.*Y./(Y+2*G));
z = sum(VF.*Y./(Y+2*G));
u = sum(VF./(Y+2*G));
v = sum(VF./G);
w = sum(VF.*G);
A = 4*x + z*z/u;
B = 2*y + z*z/u;
E = 1/u;
F = z/u;
D = 1/v;
M = w;
c = zeros(6,6);
c(1,1) = A;
c(1,2) = B;
c(1,3) = F;
c(2,1) = B;
c(2,2) = A;
c(2,3) = F;
c(3,1) = F;
c(3,2) = F;
c(3,3) = E;
c(4,4) = D;
c(5,5) = D;
c(6,6) = M;
end

[Walter Roberson]

Sig_D_P = S_A / ((kappa_f - L) * phi + (kappa_a - N));

S_A is 6 x 6

kappa_f is 1 x 3.

L is unknown size.

phi is scalar.

kappa_a is scalar.

N is unknown size.

In order for kappa_f - L to be valid, L must be scalar, or L must be 1 x 3, or L must be something by 1 or something by 3. The results of the subtraction would be 1 x 3, 1 x 3, or something by 3, or something by 3

With phi being a scalar, * phi gives the same size as the left of the * operator, so (kappa_f - L) * phi must be 1 x 3, 1 x 3, or something x 3, or something by 3

With kappa_a being a scalar, and N being an unknown size, the result of the subtraction is always valid. So (kappa_a - N) is possibly scalar, but could just as easily be 17 x 49 x 14 as far as operand consistency is concerned.

For (kappa_f - L) * phi to be consistent with (kappa_a - N), then the (kappa_a - N) part could be scalar, or could be 1 x 3, or something by 1, or something by 3, with N being scalar, or 1 x 3, or something by 1 or something by 3 respectively.

The overall expression ((kappa_f - L) * phi + (kappa_a - N)) must work out to 1 x 3 or something by 3.

There is no circumstances under which the / operator can be used between a 6 x 6 and a 1 x 3 or something by 3.

There is no circumstances under which the ./ operator can be used between a 6 x 6 and a 1 x 3 or something by 3.

Therefore S_A / ((kappa_f - L) * phi + (kappa_a - N)) and S_A ./ ((kappa_f - L) * phi + (kappa_a - N)) are both impossible operations.

[Torsten]

Maybe you want to give weights to the fitting of the A(i),B(i),C(i),D(i) and F(i) according to the number of times they appear in the matrix Sig_D(:,:,i). But in principle, I arrive at the same code as before.

Or do you want to use a specific matrix norm to compare the numel(Kf) data matrices Sig_D(:,:,i) with the simulated matrices Sig_D_P(:,:,i) ?

K1 = 60;
K2 = 40;
U1 = 30;
U2 = 20;
pt1 = 0;
pt2 = 0.4;
Kf = linspace(1, 7, 7);
Kf1 = 0;
Xl1 = 0.5;
Xl2 = 1-Xl1;
Uf1 = 0;
Uf2 = 0;
phi = Xl1 * pt1 + Xl2 * pt2;
[~,Kub_frame_1,~,Uub_frame_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf1,Uf1);
[~,Kub_frame_2,~,Uub_frame_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf1,Uf1);
Y1_Frame = Kub_frame_1 - (2*Uub_frame_1)/3;
Y2_Frame = Kub_frame_2 - (2*Uub_frame_2)/3;
[Cij_A_Frame] = Backus_average([Xl1,Xl2],[Uub_frame_1,Kub_frame_2],[Y1_Frame,Y2_Frame]);
Sij_A_Frame = inv(Cij_A_Frame);
Sig_D  = zeros(6, 6, numel(Kf));
for i = 1:numel(Kf)
    Kf2 = Kf(i);
    
    [~,Kub_ud_1,~,Uub_ud_1] = Hashin_Shtrikhman_full(K1,U1,1,pt1,Kf2,Uf2);
    [~,Kub_ud_2,~,Uub_ud_2] = Hashin_Shtrikhman_full(K2,U2,1,pt2,Kf2,Uf2);
    
    Y1_Sat = Kub_ud_1 - (2*Uub_ud_1)/3;
    Y2_Sat = Kub_ud_2 - (2*Uub_ud_2)/3;
    
    [Cij_A_Ud] = Backus_average([Xl1,Xl2],[Uub_ud_1,Uub_ud_2],[Y1_Sat,Y2_Sat]);
    Sij_star = inv(Cij_A_Ud);
    Sig_D(:,:,i) = Sij_A_Frame - Sij_star;
end
% solve this function for unknowns A, B, F, C, D, N, L
%Sig_D_P = Sig_D; 
%Sig_D_P = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L);
A = squeeze(Sig_D(1,1,:)).';
B = squeeze(Sig_D(1,2,:)).';
C = squeeze(Sig_D(3,3,:)).'; 
D = squeeze(Sig_D(4,4,:)).';
F = squeeze(Sig_D(3,1,:)).';
kappa_f = 1./Kf;
kappa_a = trace(Sij_A_Frame);
ydata = [A,B,C,D,F];
ysim = @(A1,B1,C1,D1,F1,LN) [A1./(kappa_f * phi + kappa_a - LN),...
                             B1./(kappa_f * phi + kappa_a - LN),...
                             C1./(kappa_f * phi + kappa_a - LN),...
                             D1./(kappa_f * phi + kappa_a - LN),...
                             F1./(kappa_f * phi + kappa_a - LN)];
ysim = @(p)ysim(p(1),p(2),p(3),p(4),p(5),p(6));
F = @(p) ysim(p) - ydata;
sol = lsqnonlin(F,[2 3 4 5 6 7])

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

sol = 1×6

0.0270 -0.0142 0.0683 0.2732 -0.0192 11.7901

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

function Sig_D_P = sovle_this(kappa_f,kappa_a,phi,A,B,F,C,D,N,L)
S_A = zeros(6,6);
S_A(1,1) = A;
S_A(2,2) = A;
S_A(1,2) = B;
S_A(2,1) = B;
S_A(1,3) = F;
S_A(2,3) = F;
S_A(3,1) = F;
S_A(3,2) = F;
S_A(3,3) = C;
S_A(4,4) = D;
S_A(5,5) = D;
S_A(6,6) = (1/2) * (A-B);
Sig_D_P = S_A / ((kappa_f - L) * phi + (kappa_a - N));
end

function [Klb,Kub,Ulb,Uub]=Hashin_Shtrikhman_full(K,U,Vfr,P,Kf,Uf)
Umin=min(U);
Umax=max(U);
Kmin=min(K);
Kmax=max(K);
if P==0
    Klb = 1./((1-P).*sum(Vfr./K));
elseif P>0
    Klb = 1./((P./Kf)+((1-P).* sum (Vfr./K)));
end
Kub = ((P * (1./(Kf+((4/3)*Umax))) + (1-P) * sum(Vfr./(K+((4/3).*Umax))))^(-1)) - (4/3)*Umax;
Zmax =(Umax/6)*((9*Kmax+8*Umax)/(Kmax+2*Umax));
if P ==0
    Zmin =(Umin/6)*((9*Kmin+8*Umin)/(Kmin+2*Umin));
elseif P>0
    Zmin =(Uf/6)*((9*Kf+8*Uf)/(Kf+2*Uf));
end
Uub = (( (P/Zmax) + (1-P) * sum (Vfr./(U+Zmax)))^(-1))-Zmax;
Ulb = (( (P/Zmin) + (1-P) * sum (Vfr./(U+Zmin)))^(-1))-Zmin;
end 
function [c]=Backus_average(VF,G,Y)
x = sum(VF.*G.*(Y+G)./(Y+2*G));
y = sum(VF.*G.*Y./(Y+2*G));
z = sum(VF.*Y./(Y+2*G));
u = sum(VF./(Y+2*G));
v = sum(VF./G);
w = sum(VF.*G);
A = 4*x + z*z/u;
B = 2*y + z*z/u;
E = 1/u;
F = z/u;
D = 1/v;
M = w;
c = zeros(6,6);
c(1,1) = A;
c(1,2) = B;
c(1,3) = F;
c(2,1) = B;
c(2,2) = A;
c(2,3) = F;
c(3,1) = F;
c(3,2) = F;
c(3,3) = E;
c(4,4) = D;
c(5,5) = D;
c(6,6) = M;
end